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Introduction 


In  the  application  of  composite  materials  to  structural  design,  the  accurate  prediction  of  residual 
strength  is  of  prime  interest  to  design  engineers  to  ensure  that  the  material  efficiency  afforded 
by  composites  meet  service  load  conditions  while  demonstrating  survivability  and  damage  toler¬ 
ance  to  projected  damage-causing  threats  in  the  intended  service  environment.  While  providing 
more  efficient  structures,  composites  are  inherently  more  complex  material  systems  than  metallics 
and  require  more  elaborate  mathematical  models  to  simulate  their  behavior.  The  mathematical 
representation  of  composite  structures  is  complicated  by  the  dimensional  scale  at  which  various 
phenomena  are  modeled.  Three  modeling  domains  may  be  loosely  defined  as  the  microscopic, 
lamina  and  component  level.  At  a  microscopic  level  of  analysis,  the  micromechanical  considera¬ 
tion  of  fiber  and  matrix  characteristics  are  predominant.  These  include  fiber  waviness  and  weave 
type,  defects  such  as  fiber/matrix  debonding,  fiber  breakage,  matrix  cracking,  porosity  and  inter¬ 
active  fracture  phenomena.  It  is  within  this  regime  that  material  failure  criteria  are  developed 
from  the  analysis  of  micromechanical  damage  modes.  At  the  lamina  level  of  consideration,  the 
assumption  of  a  continuum  may  be  applied  to  represent  a  composite  lamina  as  a  continuous  ho¬ 
mogeneous  material  with  associated  material  property  tensors  defining  basic  elastic  stress/strain 
relations  derived  from  en^neering  mixture  relations  of  the  bulk  properties  of  the  constituent  fiber 
and  matrix  materials.  Ply  thickness  variations,  failure  criteria,  material  damage  laws,  debonding 
and  delamination  growth  between  adjacent  laminae  and  material  nonlinearity  lead  to  mathematical 
representations  that  model  the  progression  of  damage  and  failure  of  individual  plies.  At  the  high¬ 
est  level  of  representation,  the  mathematical  modeling  of  a  composite  component  generally  must 
account  for  complex  geometries,  applied  loading  and  support  conditions.  Laminate  properties  may 
be  represented  by  a  summation  of  averaged  ply  properties  and  assumed  as  continuous  throughout 
the  local  material  domain.  Depending  on  the  geometric  complexity  of  the  problem  being  studied, 
various  approaches  may  be  used  to  model  the  problem.  Mathematical  modeling  may  utilize  exact 
formulations  based  on  elasticity  theory,  structural  engineering  approximations  such  as  beam,  plate 
and  shell  assemblages,  or  numerical  approaches  such  as  finite  difference  or  finite  element  theory. 
Regardless  of  the  modeling  approach,  the  global  response  of  the  component  to  applied  loading  is 
sought.  This  response  includes  the  basic  displacement  and  internal  stress  distribution  from  which 
local  material  failure  modes  and  ply  instabilities  may  be  predicted.  The  analytical  prediction  of 
residual  strength  ranges  over  three  basic  scale  regimes,  each  of  which  requires  different  modeling 
considerations. 

A  myriad  of  analytical  and  numerical  approaches  for  predicting  residual  strength  in  composites 
have  been  formulated  by  various  researchers.  These  analyses  have,  in  general,  been  restricted  to 
selected  damage  modes  and  to  idealized  geometric  configurations.  An  overview  of  selected  anal¬ 
yses  include  those  based  on  fracture  mechanics,  structural  engineering  approximations  and  finite 
element  theory.  Reference  [1]  presents  a  fracture  mechanics  approach  in  which  the  observed  impact 
damage  zone  is  replaced  by  an  elliptical  edge  crack.  In  this  model  it  is  assumed  that  the  crack  grows 
laterally  across  the  specimen  width  without  increasing  in  depth.  Under  increasing  load,  shearing 
at  the  base  of  the  crack  effectively  separates  or  isolates  this  layer  from  the  rest  of  the  laminate  and 
precipitates  ultimate  tensile  failure.  Similar  approaches  have  been  applied  in  which  damage  due  to 
impact  is  approximated  as  through- thickness  holes,  cracks  or  inclusions  [2-5].  Failure  of  laminates 
due  to  applied  in-plane  compressive  loading  have  been  studied  assuming  delamination  growth  and 
buckling  failure  as  the  primary  failure  modes.  Reference  [6]  considers  a  multi-beam  model  to  rep¬ 
resent  delaminations  and  assesses  the  effects  of  support  conditions,  planarity  assumptions  of  beam 
rotation  and  post-buckled  behavior  on  total  load  carrying  capability.  Reference  [7]  adopts  a  beam 
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model  to  assess  delamination  growth,  stability  and  arrest  in  a  1-D  model  of  laminate  failure.  The 
effect  of  bridging  due  to  transverse  normal  stitching  on  delamination  growth  is  presented  in  Refer¬ 
ence  [8].  The  study  of  delamination  in  a  homogeneous,  axially  loaded  plate  is  presented  in  Reference 
[9]  in  which  the  effect  of  delamination  position,  size,  thickness  and  assumed  boundary  conditions  on 
critical  buckling  loads  are  considered.  The  effect  of  an  elliptical  delamination  in  a  quasi-isotropic 
composite  plate  is  studied  in  Reference  [10].  The  delamination  is  assumed  to  be  near-surface  and 
a  Rayleigh-Ritz  approach  used  to  assess  the  effect  of  delamination  shape,  orientation  and  layup 
on  critical  buckling  strains.  A  study  of  stability  in  random,  short-fiber  composite  plates  is  pre¬ 
sented  in  References  [11,  12]  which  also  employ  a  Rayleigh-Ritz  analysis  and  assess  delamination 
growth  and  coupled  local/global  buckling  phenomena.  General  studies  on  delamination  stability 
such  as  [13,  14]  demonstrate  that  delamination  growth  can  proceed  sublaminate  buckling.  Finite 
element  analysis  of  residual  strength  in  laminates  based  on  assumed  progressive  in-plane  material 
failure  have  been  formulated  utilizing  specialized  ply  damage  laws.  Bolted  composite  joints  are 
analyzed  in  Reference  [15]  assessing  tension  and  shear-out  failure  modes.  Reference  [16]  considers 
a  pin-loaded  hole  in  relation  to  in-plane  material  failure  modes.  Laminated  composites  containing 
stress  concentration  due  to  a  circular  cut-out  is  analyzed  in  Reference  [17].  In  the  above  analyses 
of  in-plane  failure,  unique  material  damage  laws  are  developed  which  differentiate  between  fiber 
and  matrix  modes  and  includes  nonlinear  shear  stress  behavior.  The  residual  tensile  strength  of 
impacted  laminates  is  predicted  using  a  ply-level  analytical  formulation  in  [18].  The  developed 
analysis  avoids  plate  assumptions  of  plane  stress  and  utilizes  a  fully  3-D  finite  element  represen¬ 
tation  of  the  laminate.  Damage  modes  are  limited  to  delaminations  and  fiber  breakage  which  are 
represented  as  a  through-thickness  line  crack.  A  progressive  failure  analysis  is  presented  in  Refer¬ 
ence  [19]  in  which  the  layerwise  laminate  theory  of  [20]  is  used  to  model  the  laminate  continuum. 
This  analysis  improves  upon  other  progressive  failure  analyses  such  as  those  presented  in  References 
[21-24]  by  including  failure  through  delamination.  The  Tsai-Wu  failure  criteria  is  used  to  predict 
ply  failure  and  a  stiffness  reduction  model  is  used  to  approximate  post  failure  ply  behavior.  Axial 
tensile  loading  is  assumed  and,  in  comparison  to  other  structural  engineering  approximations,  it  is 
concluded  that  a  fully  3-dimensional  analysis  is  required  for  accurate  failure  prediction. 

The  purpose  of  the  present  investigation  is  to  identify  the  various  failure  phenomena  involved 
in  damaged  composite  laminates  and  to  discuss  mathematical  modeling  approaches  that  have  been 
developed  to  represent  such  phenomena.  It  must  be  emphasized  that  research  into  composite  ma¬ 
terials  has  given  rise  to  an  extensive  literature  on  composite  behavior.  The  present  effort  is  limited 
to  providing  a  representational  survey  of  various  topics  and  is  constrained  to  a  selective  set  of  ana¬ 
lytical  developments  with  which  to  broadly  outline  the  prediction  of  residual  strength  in  composites 
structures.  When  possible,  extensions  to  overcome  limitations  in  existing  analyses  are  suggested. 

The  following  sections  discuss  modeling  considerations  which  include  a  brief  description  of  nonde¬ 
structive  evaluation  techniques  (NDE)  currently  used  in  determining  initial  damage  states,  charac¬ 
teristics  of  composite  laminae  and  laminate  response  to  applied  loads.  Finally,  details  of  a  proposed 
numerical  finite  element-based  approach  are  presented  in  which  all  important  damage  modes  and 
their  interaction  may  be  accounted  for  through  appropriate  solution  algorithms.  Such  an  approach 
is  deemed  capable  of  incorporating  all  relevant  phenomena  in  a  predictive  methodology  to  accu¬ 
rately  determine  the  ultimate  strength  of  composite  materials  in  arbitrary  structural  geometries 
under  general  loading  and  support  conditions. 
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Nondestructive  Damage  Characterization 


Regardless  of  how  accurate  any  developed  analytic  methodology  is  in  representing  a  damaged  com¬ 
posite  laminate,  the  prediction  of  residual  strength  is  entirely  dependent  on  the  accuracy  of  the 
input  defining  the  internal  damage.  This  accuracy  is,  in  turn,  dependent  on  the  capability  of  ex¬ 
isting  inspection  techniques  to  resolve  internal  damage  and  distinguish  different  damage  modes. 

Inspection  procedures  may  be  destructive  or  nondestructive  in  nature.  Destructive  characterized 
tion,  in  which  a  damaged  laminate  is  physically  sectioned  and  microscopically  inspected,  is  useful 
in  providing  a  quantified  assessment  of  laminate  damage  type  and  location  due  to  specific  impact 
events,  environmental  effects  or  manufacturing  defects.  Such  an  examination,  however,  precludes 
subsequent  experimental  testing  of  the  actual  specimen  in  measuring  residual  load— carrying  capa¬ 
bility  and  a  post-failure  examination  of  failure  modes.  Nondestructive  evaluation  (NDE)  techniques 
provide  detailed  information  of  initial  defects  and  internal  damage  while  allowing  subsequent  ex¬ 
perimental  assessment  of  residual  strength.  This  information  used  in  mathematical  simulation 
of  residual  strength  make  possible  accept/reject  decisions  during  fabrication  or  the  assessment  of 
repair/no-repair  options  for  composite  components  effected  in  service  applications. 

Several  survey  articles  of  current  state-of-the-art  NDE  techniques  may  be  found  in  References 
[25-28].  The  aim  of  these  techniques  is  to  obtain  information  or  ‘signatures’  of  internal  structural 
features  which  are  used  to  locate  and  differentiate  various  defects  and  damage.  Such  internal 
features  include  manufacturing  defects  such  as  voids,  porosity,  cracks  and  ply  debonding,  service 
environmental  effects  such  as  corrosion  and  moisture  levels,  and  damage  due  to  overloads  or  impacts. 

Of  the  multitude  of  NDE  techniques  available  as  mature  technologies,  a  brief  outline  of  several 
versatile  methods  are  outlined  below  to  give  an  indication  of  current  capabilities  in  characterizing 
internal  damage.  The  selected  techniques  which  are  currently  in  wide  use  are  acoustics/ultrasonics, 
thermography  and  radiography. 


Acoustics/Ultrasonics 


Within  this  class  of  techniques,  sound  or  pressure  pulses  are  used  to  obtain  information  on  internal 
damage  [29-35].  One  technique  utilizes  the  acoustical  emissions  recorded  during  experimental 
testing  at  the  onset  of  component  failure  to  correlate  sound  frequencies  with  specific  failure  modes 
such  as  delaminations,  matrix  cracking  and  fiber  breakage.  The  simplest  use  of  acoustical  emissions 
is  the  Tap  Test  in  which  a  rigid  object  is  lightly  tapped  along  the  outer  surface.  The  change  in  sound 
pitch  of  the  vibrating  laminate  can  indicate  the  presence  of  near  surface  damage,  predominantly 
delaminations.  Ultrasonic  inspection  techniques  utilize  a  transducer  to  introduce  pressure  pulses 
into  the  material.  The  governing  equations  used  to  model  waves  in  an  isotropic  solid  are  given  by 
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where  u  and  v  are  the  displacements  in  the  x  and  y  directions,  t  is  time,  A  and  fi  are  the  Lame 
constants,  and  p  is  the  density.  Internal  damage  will  interact  and  scatter  the  stress  waves  which 
are  recorded  as  backscatter  or  through-transmission  changes  in  the  input  signal.  This  technique 
can  detect  defects  and  damage  modes  such  as  voids  and  regions  of  material  porosity  in  addition  to 
delaminations,  regions  of  matrix  cracking  and  fiber  breakage.  The  resolution  of  damage  diminishes 
with  increasing  specimen  thickness,  attenuating  features  inherent  to  the  material  and  an  increasing 
number  of  damage  sites  along  the  direction  of  the  wavefront.  An  example  of  ultrasonic  imaging 
is  presented  in  Figure  1  showing  the  profile  of  square  metallic  shim  positioned  at  the^  midplane 
of  a  1”  thick  S-2  Epoxy  woven  laminate.  The  figure  demonstrates  the  capability  of  imaging  a 
deeply  imbedded  defect  in  a  woven  cloth  laminate  which  presents  significant  attenuation  of  induced 
pressure  pulses  due  to  the  complex  fiber  arrangement  in  the  individual  ply  weaves. 


Figure  1.  Ultrasonic  image  of  a  thick  composite  panel  with  simulated  damage. 


4 


Thermography 


Thermal  imaging  of  internal  damage  is  based  on  the  change  of  temperature  distribution  in  a  mate¬ 
rial.  A  common  form  of  thermography  utilizes  the  change  in  thermal  conductivity  due  to  material 
discontinuities  associated  with  different  damage  states  [36,37].  The  governing  equation  may  be 
given  by 


6T  8  8T.  ,  8  ,,  8T^  8  8T , 

“'’ft  = 


(2) 


where  T  is  the  temperature,  p,  c  and  k  are  the  density,  specific  heat  and  thermal  conductivity, 
respectively.  In  applying  the  technique,  thermal  energy  is  introduced  into  the  specimen  and  in¬ 
frared  thermal  emissions  are  measured  to  construct  images  of  regions  of  nonuniform  heat  flow  thus 
profiling  the  damage.  Thermal  gradients  may  be  measured  on  the  same  or  opposite  side  of  applied 
heating.  Image  resolution  declines  rapidly  with  increasing  depth  of  internal  damage.  Alternatively, 
thermograms  may  be  generated  by  imaging  the  internal  temperature  changes  caused  by  applied 
loading.  A  fundamental  relationship  exists  relating  internal  stresses  to  generated  temperature  dis¬ 
tributions  which  are  quantified  through  the  equation  of  state  of  the  material  being  studied.  This 
phenomenon  is  formalized  in  the  mathematical  theory  of  thermoelasticity.  For  an  isotropic  mate¬ 
rial,  the  thermoelastic  equilibrium  equations  relating  stresses  and  temperature  are  given  by 


8x  8y  8z 

^’’^xy  ^^yy  ,  ^'^yz 

8x  8y  8z 

8Txz  ^"^yz  .  zz 

8x  8y  8z 


aE  8T 
l-2iy8x 
aE  8T 
I  —  2i/  8y 
aE  8T 
I  —  2iy  8z 


(3) 


where  E  is  the  Young  s  modulus,  u  is  the  Poissons  ratio  and  a  is  the  coefficient  of  thermal  expansion. 
Current  sensors  can  detect  temperature  changes  on  the  order  of  O.OOrC  which  gives  a  resolution 
comparable  to  common  strain  gauges.  As  stated  in  Reference  [25],  refinements  in  thermographic 
imaging  may  be  expected  through  the  development  of  more  sensitive  signal  analysis  tools.  An 
example  of  thermographic  imaging  is  presented  in  Figure  2  showing  the  combined  stress  profile 
around  a  1/4”  diameter  open  hole  in  a  unidirectional  Aramid  fiber  laminate  under  uniform  applied 
tensile  loading. 
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Figure  2.  Thermograph  of  stresses  around  an  open  hole  m  a  loaded  specimen. 


Radiography 

Radiography  is  a  powerful  technique  for  characterizing  internal  damage  in  composites  [38,  39,  40]. 
A  common  form  of  this  technique  utilizes  high  energy  electromagnetic  radiation  in  the  form  of  x- 
rays  to  measure  the  changes  in  absorption  due  to  variations  in  density.  Radiography  is  an  extremely 
sensitive  technique  and  yields  high  resolution  images  of  internal  damage.  Current  reso  utions  ob¬ 
tainable  are  on  the  order  of  0.1%  changes  in  density.  The  resulting  2-D  damage  Vromes^y  he 
further  refined  using  various  image  enhancement  techniques.  In  Computed  Tomography  (CT),  l-D 
radiographic  measurements  are  made  along  rays  through  the  specimen  at  a  large  number  of  points 
around  the  perimeter.  Mathematical  image  reconstruction  algorithms  are  then  utihzed  to  genera  e 
a  2-D  image  of  damage  in  the  plane.  This  may  be  repeated  through  the  thickness  of  the  specimen  to 
obtain  a  fuUy  3-D  image  of  internal  damage.  This  capability  coupled  with  the  high  sensitivity  and 
resolution  of  radiography  provides  a  valuable  tool  for  determining  the  3-D  distribution  of  various 

composite  damage  modes. 

Using  basic  NDE  technologies  such  as  those  mentioned  above  to  obtain  information  of  internal 
damage,  various  image  processing  algorithms  are  available  for  extracting  information  and  charac¬ 
terizing  damage  [41-44].  Reference  [45]  explores  the  post-processing  of  reconstructed  images  from 
computed  tomography  to  identify  the  3-D  distribution  of  internal  delaminations  fiber  and  matnx 
damage  regions  through  modal  separation.  Images  were  obtained  from  thick  S-2  Glass  composite 
panel? ballistically  impacted  in  a  study  presented  in  Reference  [46].  Figure  3  shows  a  representati 
CT  image  of  impact  induced  damage  at  a  particular  depth  level.  Research  has  focused  on  a 
matically  interpreting  the  damage  signature  present  in  the  i^mage  to  differentiate  various  dam^e 
modes.  Figure  4  shows  a  histogram  of  pixel  intensities  which  exhibits  various  clusters  of  nor  y 
distributed  intensity  levels.  These  local  distributions  or  modal  patterns  are  logically  conjectured 
to  represent  separate  damage  modes  and  automated  algorithms  are  being  developed  to  interpret 

these  damage  signatures. 
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Figure  3.  CT  image  of  damage  in  a  composite  laminate. 


Figure  4.  Histogram  of  image  intensities  showing  damage  modes. 


The  objective  of  the  image  analysis  is  to  extract  maximum  information  about  the  internal  damage 
state  which  is  of  paramount  importance  as  input  to  any  developed  residual  strength  prediction 
methodology. 


Characteristics  of  Composite  Laminae 


In  the  mathematical  representation  of  a  composite  component,  the  primary  constituents  of  a  lami¬ 
nate  are  the  plies  or  laminae.  A  lamina  consists  of  an  aggregate  of  high  strength  fibers  embedded  in 
a  surrounding  polymeric  organic  matrix.  Plies  may  be  formed  with  long  continuous  fibers  aligned 
in  a  principle  direction  to  constitute  a  tape,  with  short,  randomly  oriented  fibers  forming  a  trans¬ 
versely  isotropic  layer,  or  with  fibers  arranged  in  a  multitude  of  different  weave  patterns  constituting 
a  woven  cloth  ply.  This  biphasic  material  system  is  characterized  by  fiber  orientation,  distribution 
and  material  properties,  matrix  material  properties  and  numerous  micromechanical  aspects  includ¬ 
ing  initial  manufacturing  defects  and  service  load  induced  damage. 

In  the  development  of  a  residual  strength  prediction  methodology,  several  aspects  of  composite 
laminae  are  of  critical  importance  in  the  development  of  accurate  models  which  include  the  mathe¬ 
matical  representation  of  a  ply  as  a  homogeneous  continuum,  nonlinear  material  properties,  failure 
theories,  ply  strengths  and  damage  laws.  These  issues  are  discussed  in  the  following  sections. 


Continuum  Representation 

A  fundamental  aspect  in  the  mathematical  representation  of  laminated  composites  is  the  treatment 
of  the  ply.  The  density  and  distribution  of  fibers  in  a  ply  is  altered  during  curing  by  the  flow  of 
resin  and  may  be  described  statistically.  Fiber  volume  and  the  nature  of  the  fiber  bonding  with 
the  surrounding  matrix  is  important  in  approximating  the  gross  behavior  of  the  laminae  [47].  The 
modeling  of  plies  or  laminae  involve  micromechanical  considerations  of  fiber  continuity,  waviness, 
bonding  characteristics  with  the  surrounding  matrix,  the  effect  of  broken  fibers  on  local  stress 
redistribution  and  fracture  phenomenon  dealing  with  the  propagation  and  interaction  of  fracture 
surfaces.  For  the  purpose  of  macroscopic  analysis,  these  effects  are  modeled  using  a  continuum 
assumption  in  which  the  lamina  elastic  properties  are  determined  from  a  simple  mixture  rule  and 
the  bulk  material  properties  of  the  fiber  and  matrix.  Micromechanical  defects  such  as  fiber-matrix 
debonding,  fiber  breakage,  matrix  cracking,  voids,  porosity,  moisture  and  temperature  effects  are 
represented  through  appropriate  changes  to  the  primary  variables  -  moduli  -  describing  the  elastic 
continuum.  These  micromechanical  bases  for  macroscopic  behavior  may  be  determined  analytically 
by  assuming  a  representative  unit  cell  or  volume  element  which  consists  of  an  idealized  arrange¬ 
ment  of  fibers  in  a  surrounding  matrix  and  incorporates  assumed  micromechanical  damage  modes. 
This  unit  cell  is  assumed  as  a  fundamental  description  of  the  ply  continuum  and  elasticity  theory 
is  used  to  determine  the  elastic  behavior  of  this  fundamental  ply  unit.  An  example  of  a  typical 
two-dimensional  unit  cell  is  depicted  in  Figure  5  showing  idealized  distributions  of  broken  fibers 
with  regions  of  fiber /matrix  debonding.  The  results  of  the  stress  and  failure  analysis  of  the  unit 
cell  are  scaled  up  to  approximate  macroscopic  ply  behavior. 
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Figure  5.  Typical  unit  cell  representation  of  ply  continuum  in  micromechanical  analysis. 

The  idealization  of  a  ply  as  an  elastic  continuum  inevitably  contributes  to  the  disparity  between 
analytical  predictions  and  actual  laminate  behavior.  The  basis  for  this  disparity  is  inherent  to  the 
assumption  of  homogeneity  in  a  typical  ply  and  may  be  contrasted  with  the  actual  discontinuous 
and  heterogenous  ply  composition  as  depicted  in  Figure  6  which  shows  a  photomicrograph  from 
Reference  [48]  of  a  ply  with  a  central  vertical  crack. 


Figure  6.  Photomicrograph  of  a  typical  ply  with  a  vertical  crack. 
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The  complexity  of  the  underlying  material  composition  of  individual  plies  makes  the  mathematical 
treatment  approximate  or  statistical  in  nature.  Thus,  significant  modeling  assumptions  are  present 
in  the  basic  representation  of  composite  laminae  which  introduce  an  inaccuracy  and  variability  to 
any  prediction  of  composite  behavior. 


Nonlinear  Ply  Properties 


Few  materials  maintain  a  linear  elastic  or  Hookian  response  to  failure,  yet  the  nonlinearities  exhib¬ 
ited  near  the  point  of  failure  for  many  materials  are  negli^ble.  For  composite  materials,  nonlinear 
viscoelastic  and  plastic  behavior  is  commonly  neglected  and  most  elastic  components  may  be  repre¬ 
sented  as  Hookian  or  exhibiting  a  weak  nonlinear  elastic  reponse.  Typical  fibers  used  in  composite 
material  systems  are  high  strength  yet  brittle  and  exhibit  a  linear  stress-strain  relation.  Matrix 
materials  are  usually  compliant  and  exhibit  various  degress  of  nonlinear  response  in  high  stress 
regimes.  As  depicted  in  Figure  7,  the  combined  effect  on  lamina  behavior  is  commonly  linear  with 
increasing  stress  until  fiber  failures  begin.  In  addition,  lamina  failure  strain  is  typically  close  to  the 
fiber  failure  strain. 


STRAIN 

Figure  7.  Stress-strain  behavior  of  ply  constituents. 
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For  weakly  nonlinear  ply  elastic  components,  an  approximation  to  the  nonlinearities  may  be  made 
through  a  bilinear  stress-strain  relation  as  shown  in  Figure  8.  The  use  of  a  bilinear  relationship 
greatly  simplifies  analytic  formulations  over  higher-order  functional  representations. 


a 


Figure  8.  Bilinear  stress-strain  law. 


Certain  composite  material  systems  can  exhibit  significant  nonlinear  behavior  in  all  or  specific 
elastic  moduli.  Experimental  verification  of  unidirectional  nonlinear  ply  response  under  tensile 
loading  is  presented  in  Reference  [49].  A  gross  measure  of  nonlinearity  is  presented  in  the  form  of 
an  empirical  relationship  between  an  averaged  Young’s  modulus  and  stress  level  given  by 


E  =  Eo  +  21(7  (4) 

For  the  laminate  considered,  the  increase  in  Young’s  modulus  from  the  initial  zero  stress  state 
to  ultimate  tensile  strength  was  on  the  order  of  30%.  Composite  laminates  exhibit  a  significant 
nonlinearity  primarily  in  the  elastic  response  in  the  shear  modulus  [90].  As  detailed  in  Reference 
[51],  a  piece- wise  linear  approximation  is  given  by 


{Gl2)n  = 


(Ari2)n 

(A7l2)n 


(5) 


where  at  the  segment  along  the  shear  stress/strain  curve,  the  relationship  between  Ari2  and 
A712  is  given  as  a  function  of  material  compliance  and  reduced  stiffness  terms. 

A  representation  developed  in  [90,  81]  for  the  shear  modulus  is  given  by 


712  =  (7^)  +  a(ri2f 
Gi2 


(6) 
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where  Gu  is  the  inital  shear  modulus  and  a  is  an  experimentally  determined  constant.  The  tangent 
shear  modulus  at  any  position  on  the  shear  stress-strain  curve  is  given  by 

712  =  =  (^)  +  3a(ri2)^  (7) 

It  has  been  noted  in  [17,  16,  51,  52]  that  the  nonlinearity  does  not  become  significant  until  the  ply 
shear  stresses  become  comparable  to  the  lon^tudinal  tensile  stresses.  In  an  analytical  methodol¬ 
ogy,  the  nonlinear  expression  may  be  used  in  an  iterative  solution  algorithm  for  predicting  laminate 
response  with  increasing  load.  To  avoid  the  computational  expense  of  a  nonlinear  solution  scheme, 
the  shear  modulus  may  be  taken  as  the  initial  modulus,  a  median  tangent  modulus,  or,  frequently, 
as  a  linear  approximation  using  a  secant  modulus  defined  by  the  shear  strength  and  ultimate  shear 
strain.  The  various  definitions  of  the  shear  modulus  are  depicted  in  Figure  9. 


Depending  on  the  laminate  ply  layup  and  magnitude  of  the  shear  stress,  the  various  approxi¬ 
mations  in  linear  analysis  may  be  unacceptably  inaccurate  in  predicting  through-thickness  stress 
distribution  and,  hence,  in  predicting  ply  level  damage. 


Ply  Failure  Criteria 


A  primary  consideration  in  the  analysis  of  composite  structures  is  the  accurate  prediction  of  ply 
failure  due  to  a  specific  state  of  stress.  The  complexity  of  the  biphasic  ply  material  makes  the 
analysis  of  local  failure  difficult.  Under  applied  external  loading,  the  local  stress  flow  yields  a  large 
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number  of  identifiable  failure  interactions  and  damage  modes,  any  one  of  which  can  become  a  crit¬ 
ical  ‘weak-link’  in  the  prediction  of  ply  failure.  The  difficulty  in  this  endeavor  is  aptly  expressed 
by  a  passage  from  [53]:  “The  microstructural  aspects  of  failure  are  of  such  complexity  that  there 
is  little  hope  of  resolution  of  this  problem  on  the  basis  of  micromechanics  methods.  Such  meth¬ 
ods  would  require  analytical  detection  of  successive  microfailures  in  terms  of  microstress  analysis 
and  microfailure  criteria  and  prediction  of  the  coalescence  of  some  of  them  to  form  macrofailures 
which  is  an  intractable  task.”  To  make  the  mathematical  modeling  of  such  phenomena  tractable, 
simplifying  assumptions  are  made  utilizing  a  fracture  mechanics  or  continuum  approach.  In  a 
fracture  mechanics  framework,  failure  of  a  ply  is  based  on  the  predicted  growth  of  initial  cracks. 
A  continuum  approach  assumes  that,  mathematically,  accurate  failure  criteria  involving  homoge¬ 
neous  elastic  constants  exist  at  a  macroscopic  level.  Material  failure  mechanisms  have  been  well 
established  for  isotropic  materials  in  which  different  definitions  for  failure  are  used  such  as  ultimate 
strength,  yielding,  proportional  limit,  endurance  limit  and  maximum  working  stress.  Numerous 
failure  theories  have  been  advanced  by  Rankine,  Tresca-Guest,  St.  Venant,  Beltrami,  Huber-Von 
Mises-Hencky,  and  Mohr  [54].  These  theories  are  regarded  as  single-parameter  criteria  as  they  are 
based  on  a  single  measure  of  failure,  usually  the  yield  stress  of  the  material.  These  criteria  are, 
however,  inadequate  to  represent  the  failure  of  more  complex  materials  in  which  strength  differs 
in  tension  and  compression  and  stress  interaction  effects  are  significant.  Subsequent  failure  criteria 
have  been  developed  for  anisotropic  materials  from  simple  criteria  such  as  maximum  stress/ strain 
to  those  of  Hill  [55],  Tsai-Hill  [56],  Hoffman  [57]  and  Tsai-Wu  [58]  which  are  referred  to  as  quadratic 
polynomial  theories  [59,  60,  61].  Each  may  be  deduced  from  the  general  tensor  polynomial  criteria 
given  by 

FI  =  FiOi  -t-  FijOiOj  -f  FijkCTiajak  +  -  (8) 

where  tr,-  are  stress  tensor  components  in  principle  material  directions  and  F,-,  Fij,  Fijk  are  compo¬ 
nents  of  the  strength  tensors  and  FI  is  the  failure  index.  Failure  is  predicted  when  FI  >  1.  For 
practical  applications,  the  terms  in  equation  (8)  higher  than  quadratic  are  usually  set  to  zero  due 
to  the  diminishing  returns  of  including  higher-order  polynomial  terms  in  predicting  failure  together 
with  the  increased  cost  of  experimentally  determining  these  higher-order  strength  components. 
In  2-D  problems,  a  problem  in  tensor  polynomial  criteria  is  the  ‘interaction  term’,  F^,  which  is 
constrained  by  the  following  ‘stability  criterion’  given  in  [62]  as 

F11F22  -  >  0  (9) 


It  is  argued  in  [62]  that  the  difficulty  in  experimentally  determining  the  interaction  term,  F12,  is 
highly  sensitive  to  experimental  measurement  in  regajds  to  satisfying  the  above  stability  criterion. 
It  is  demonstrated  that  the  term  may  be  dropped  without  adversely  effecting  failure  prediction.  For 
the  2-D  case,  an  approach  for  approximating  the  interaction  term  tis  a  function  of  normal  strength 
measures  is  presented  in  [63]  which  yields 

- - 1 - 

A  similar  approach  is  adopted  to  determine  the  interaction  terms  arising  in  3-D  problems. 


Of  the  various  anisotropic  failure  criteria,  the  maximum  stress  criterion  is  independent  of  the  par¬ 
ticular  failure  mode  and  ultimate  ply  failure  is  predicted  when  any  one  of  the  following  conditions 
are  satisfied 
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where  a\,  a-i  and  a-^  are  the  normal  stress  components  and  <74,  CT5  and  are  the  shear  stress 
components.  The  strength  measures  are  given  by  Xt,  Xc,  It,  Yc,  Zt  and  Zc  which  are  the  normal 
tensile  and  compressive  strengths  in  the  x,  y  and  z  directions,  respectively,  and  R,  S  and  T  are 
the  shear  strengths  in  the  yz,  xz  and  xy  planes,  respectively.  The  maximum  strain  criterion  is  of 
identical  form  as  the  above  with  stresses  replaced  by  corresponding  strains  and  strength  measures 
replaced  by  associated  maximum  strain  measures. 


The  general  Hoffman  criterion  is  given  by 

FI  =  Fiai  +  F^cTi  +  F^cr^  +  +  Fi2(Ti(T2  +  Fi^cticts  +  i^22<^2  d" 
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For  2-D  plane  problems  the  Hoffman  criterion  is  obtained  by  dropping  appropriate  terms  to  yield 
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The  general  Hill  criterion  is  based  on  exclusively  quadratic  terms  given  by 

FI  =  Fii^f  +  Fi20'i(T2  +  +  F22<^2  +  ^23<^20’3  +  Fzz(T^  +  F44<^4  +  F55<t|  +  Fg(t\  (15) 
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in  which  the  value  for  the  normal  tensile  strengths,  X,  Y  and  Z,  are  taken  as  compressive  or  tensile 
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depending  on  the  sign  of  the  corresponding  normal  stresses. 
The  general  Tsai-Wu  criterion  is  given  by 
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For  2-D  plane  problems  the  Tsia-Wu  criterion  is  presented  in  Reference  [64]  as 


FI  =  FiCTi  +  F2(T2  +  FiiCTi  +  2Fi2<Ti<T2  +  F22O2  +  ■^66<’‘i2 
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The  above  criteria  have  been  successfuly  used  to  predict  failure  in  generally  anisotropic  materials. 
One  drawback,  however,  is  that  without  modification  these  criteria  do  not  predict  the  underlying 
failure  mode  which  is  important  for  design  purposes.  In  Reference  [53]  the  basic  form  of  anisotropic 
failure  theories  are  used  to  develop  criteria  to  account  for  specific  failure  modes  of  the  basic  ply 
constituents  such  as  fiber  rupture  in  tension  or  buckling  in  compression  and  general  matrix  failure. 
Recently,  ply  failure  criteria  have  been  proposed  accounting  for  ply  nonlinearities  and  incorporated 


15 


into  successful  analyses  of  composite  failure  such  as  those  presented  in  References  [65,  66].  A  tensile 
matrix  failure  criterion  is  stated  therein  as 


^xy  _  g2 


(21) 


where  <Ty  and  (Txy  are  the  layer  transverse  tensile  stress  and  shear  stress,  respectively,  and  Yt  and  Sc 
are  the  transverse  tensile  strength  and  the  modified  in  situ  shear  strength,  respectively.  Compressive 
matrix  failure  may  be  predicted  using  the  Hashin  Failure  Criterion  [53]  which  is  stated  as 
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A  modification  of  the  Yamada-Sun  failure  criteria  presented  in  [65,  66]  may  be  used  to  predict  both 
fiber-matrix  shearing  and  fiber  breakage  as  given  by 
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where  Cx  and  Xt  are  the  longitudinal  tensile  stress  and  strength  in  each  ply,  respectively.  In  using 
the  above  criteria,  matrix  or  fiber  failure  is  predicted  when  Cmt,  ^mc  or  ey  >  1. 


Additional  failure  modes  have  been  identified  and  examined  in  the  literature.  One  such  mode 
is  a  microbuckling  phenomenon  giving  rise  to  propagating  kink  bands  in  ply  fiber  arrays  eventually 
leading  to  fiber  fracture.  As  discussed  in  References  [67,  68],  such  kink  bands  are  precipitated 
under  compressive  loading  in  primary  load-carrying  plies  usually  at  the  free  edges  which  then  tend 
to  propagate  through  the  width.  The  change  in  strain  energy,  Af/^,  of  a  unit  width  kink  band  is 
given  by 

ACf  =  {2Ajgs  -b  idsayrr.ll9,)P^]  (24) 

in  which  Aj,  dj  are  the  fiber  cross-sectional  area  and  diameter,  gj  is  the  critical  fiber  energy  re¬ 
lease  rate  corresponding  to  the  flexural  toughness  of  the  fiber,  (7ym  is  the  yield  stress  of  the  matrix 
material,  h  and  Ok  are  the  kink  zone  length  and  angle  of  rotation  with  respect  to  the  fiber  axis, 
tk  is  the  thickness  of  the  kinked  ply  layer  and  c/  is  the  fiber  volume  fraction  of  the  kinked  ply. 
Future  work  is  suggested  to  develop  a  ply  or  laminate  failure  criterion  together  with  a  damage  law 
accounting  for  kink  zone  development  and  propogation.  This  type  of  damage  is  dependent  on  the 
laminate  geometric  configuration,  loading  and  support  conditions. 

The  selection  of  appropriate  ply  failure  criteria  is  thus  an  important  aspect  of  the  mathemati¬ 
cal  modeling  of  residual  strength.  Any  lamina-level  criterion  which  is  based  on  the  assumption  of  a 
homogeneous  ply  continuum  inevitably  introduces  an  approximation  and  source  of  predictive  error 
in  the  failure  analysis  of  composite  laminates. 


Ply  Strength  Measures 

The  prediction  of  ply  failure  in  a  composite  laminate  requires  an  accurate  calculation  of  stress 
distribution  through  the  laminate  thickness  together  with  experimentally  determined  ply  strength 
measures  used  in  various  failure  criteria.  For  isotropic  materials,  a  single  strength  parameter  may 
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be  obtained  from  a  tensile,  compressive  or  shear  stress  test.  Such  materials  may  be  referred  to 
as  one-dimensional  or  one-constant  material  systems.  In  contrast,  composite  laminae  represented 
as  anisotropic  homogeneous  materials  require  several  parameters  to  characterize  the  failure  surface 
when  plotted  in  stress  or  strain  space.  These  include  compressive  and  tensile  strength  along  the 
longitudinal  and  transverse  directions  relative  to  ply  fibers,  and  shear  strengths  in  all  principle  ply 
directions.  A  typical  2-D  failure  envelope  for  a  composite  ply  is  depicted  in  Figure  10. 


Figure  10.  Typical  2-D  failure  envelope  for  a  composite  ply. 


Ply  strengths  are  dependent  on  the  basic  fiber  and  matrix  properties,  fiber  distribution  and  ori¬ 
entation,  and  on  the  integrity  of  the  interfacial  bonding  between  the  fibers  and  the  matrix  [47]. 
Under  a  general  state  of  local  stress,  a  number  of  different,  potential  weahest-links  exist  to  dic¬ 
tate  ply  strength.  Lamina  strengths  are  most  commonly  obtained  experimentally  yet  analytical 
models  may  be  developed  to  predict  strength  based  on  assumptions  regarding  the  applied  loading 
and  micromechanical  characteristics  of  the  ply  constituents.  Reference  [69]  discusses  in  detail  the 
importance  of  fiber-matrix  bonding  characteristics  and  examines  the  influence  of  bond  strength 
at  the  fiber-matrix  interface  on  ultimate  tensile  strength  by  analytically  computing  strength  mea¬ 
sures  based  on  progressive  fracturing  of  fibers.  A  micromechanical  equilibrium  element  model  is 
developed  in  which  an  initially  assumed  cylindrical  core  of  damaged  composite  is  surrounded  by  a 
concentric  array  of  undamaged  fibers.  The  result  of  the  analysis  yields  a  relationship  for  strength 
of  the  following  form 


The  derivation  of  equations  (25,26)  are  fully  detailed  and  all  variables  defined  in  Reference  [69]; 
the  equations  are  presented  to  show  the  functional  form  and  iterative  nature  of  an  analytical 
determination  of  ply  strengths  based  on  micromechanical  elasticity  models.  In  the  above,  ra,-  is  the 


17 


number  of  neaxest-neighbor  fibers  with  i  denoting  the  number  of  broken  fibers.  Qi  is  the  number 
of  broken  fiber  clusters  and  A,-  is  related  to  the  stress  concentration  in  neighboring  libers.  To 
determine  the  ply  strength,  the  applied  load  o-,-  is  incremented  until 


ffi  >  <T,+i  (27) 

which  indicates  that  the  ply  tensile  strength  is  ^ven  by  <7,-  since  additional  applied  loading  causes 
an  unstable  and  catastrophic  sequence  of  fiber  fracture  to  occur.  Other  idealized  micromechanical 
models  may  be  developed  to  analytically  estimate  other  ply  strength  measures. 

In  the  experimental  determination  of  ply  strengths,  various  tests  are  performed  using  unidirec¬ 
tional  laminates  in  which  it  is  conventionally  assumed  that  the  computed  strengths  are  fundamen¬ 
tal  constants  of  a  particular  composite  material  system.  Linear  criteria  such  as  maximum  stress 
or  maximum  strain  failure  theories  assume  an  independence  of  longitudinal,  transverse  and  shear 
strengths  thereby  simplifying  the  experimental  tests  required  to  determine  the  composite  failure 
strengths.  Quadratic  or  general  tensor  polynomial  failure  criteria  incorporate  interaction  terms 
which  require  difficult  experimental  tests  with  biaxially  loaded  test  specimens  to  determine  their 
values.  In  many  cases,  however,  it  has  been  noted  that  these  interaction  terms  may  be  neglected 
without  adversely  effecting  failure  prediction  [62]. 


An  important  observation  noted  in  Reference  [70]  states  that  the  past  assumption  of  treating 
ply  strengths  as  material  constants  is  incorrect  but  is  instead  dependent  on  laminate  thickness  and 
ply  orientation.  A  disparity  as  great  as  50%  has  been  observed  in  calculated  ply  strengths  between 
experimental  results  for  unidirectional  and  multi-directional  laminates.  Other  experimental  inves¬ 
tigations  have  corroborated  laminate  dependence  on  measured  strength  values.  In  Reference  [71] 
fiberglass  composites  demonstrated  shear  strength  values  strongly  dependent  on  laminate  thick¬ 
ness.  In  Reference  [72]  the  shear  strength  of  T300/1034-C  Graphite/Epoxy  prepregs  were  found 
to  be  a  function  of  the  volume  fraction  of  0  degree  fibers.  Further  studies  on  the  difference  of  in 
situ  strength  of  individual  plies  which  differ  from  strength  determinations  made  from  experimental 
testing  of  unidirectional  laminates  may  be  found  in  [65,  73,  74,  53,  75].  This  strength  variation  is 
attributable  to  ply  orientation,  laminate  thickness  and  residual  thermal  stresses  of  which  a  precise 
functional  relationship  is  currently  unavailable.  Empirical  relationships  are  presented  in  [65,  66] 
which  gives  expressions  for  ply  transverse  tensile  and  shear  strength  as 


Ft  =  y°(i  +  Qi 

St  =  -|-  02 


sin{A0) 

sin{A0) 


) 


) 


(28) 

(29) 


where  Yt  and  St  are  the  in  situ  transverse  tensile  and  shear  strengths.  Y°  and  S°  are  the  respective 
ply  strengths  obtained  from  the  testing  of  unidirectional  laminates.  Ad  is  the  angular  difference 
between  adjacent  ply  orientations,  N  is  the  number  of  plies  in  the  laminate  or  sublaminate  be¬ 
ing  considered  and  aj,  aj,  P\  and  P2  are  experimentally  determined  constants.  Ply  compressive 
strengths  are  assumed  to  be  of  similar  form. 


The  accurate  determination  of  ply  strengths  is  of  primary  importance  in  predicting  residual  strength 
and  further  improvements  in  testing  methods  and  analytical  models  are  needed  to  better  predict 
in  situ  ply  strength  values. 
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Ply  Damage  Models 


Once  ply  failure  is  predicted,  various  damage  laws  have  been  developed  to  account  for  degrada¬ 
tion  in  the  local  3-D  material  properties.  The  simplest  damage  law  applied  to  composite  laminae 
reduces  all  ply  properties  to  zero  if  ply  failure  has  been  predicted,  regardless  of  predicted  damage 
mode.  This  assumption  is  regarded  as  extreme  and  more  specialized  property  reductions  models 
have  been  advanced  for  the  representation  of  ply  damage. 


For  compressive  or  tensile  matrix  failure.  References  [15,  65]  have  developed  a  model  in  which 
the  in-plane  transverse  modulus,  E2  and  Poisson’s  ratios  fXij  in  the  damaged  ply  is  set  to  zero  while 
all  other  elastic  constants  are  unchanged.  This  results  in  the  ply  stiffness  matrix  given  by: 
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(30) 


Reference  [76]  discusses  a  refinement  to  the  above  model  which  may  be  regarded  as  a  limiting 
case  for  tensile  crack  density  saturation.  A  continuum  mechanics  approach  is  utilized  to  define 
material  stiffness  reductions  based  on  accumulated  transverse  crack  density.  Good  correlation  with 
experimental  data  is  obtained  although  some  minor  alterations  are  required  to  maintain  symmetry 
in  the  stiffness  matrix  [77,  78]. 


For  fiber  failure.  Reference  [79]  presents  a  model  in  which  it  is  assumed  that  in  the  damaged 
ply  the  in-plane  transverse  modulus,  E2  and  Poisson’s  ratios  /ijj  are  zero  while  the  longitudinal 
modulus,  El,  and  the  shear  modulus,  G12,  are  reduced  according  to  a  micromechanical  fiber  bundle 
failure  model  [17,  80,  81]  yielding  a  WeibuU  distribution  of  these  properties  given  by 


(31) 

1 = 

(32) 

where  E]  and  G\^  are  the  reduced  values  for  the  tensile  and  shear  moduli!,  respectively,  A  is 
the  area  of  fiber  failure  predicted  by  appropriate  fiber  failure  criteria  and  the  reference  area,  Ao,  is 
the  failure  interaction  zone  based  on  the  measured  ply  tensile  strength  Xt  [17,  80,  81,  82].  /?  is  the 
WeibuU  distribution  shape  parameter  dependent  on  the  extent  of  observed  fiber  damage.  In  the 
determination  of  the  shape  parameter,  it  has  been  noted  that  the  prediction  of  ultimate  strength  is 
somewhat  insensitive  to  the  value  of  /?  as  a  variation  of  more  than  30%  had  a  negligible  effect.  An 
alternative  approach  which  avoids  the  statistical  measurement  of  the  fiber  strength  is  presented  in 
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Reference  [19].  In  this  approach  a  stiffness  reduction  coefficient,  a,  is  utilized  to  degrade  material 
properties.  This  coefficient  is  experimentally  determined  and  is  assumed  as  a  material  property 
constant  for  the  composite  material  system  being  used.  An  ‘interactive’  degradation  model  is 
presented  based  on  a  particular  use  of  the  Tsai-Wu  failure  criterion  which  is  expressed  as 

FiOi  +  FijOiOj  =  FI  (33) 

where  FI  is  the  failure  index.  The  stiffness  reduction  due  to  damage  is  based  on  computing  the 
failure  index  for  each  component  of  stress  in  turn  to  determine  the  mode  of  damage.  Fiber  breakage 
is  assumed  if  Tnax{FI{<Ti)}  is  due  to  a\.  The  damage  law  for  this  mode  is  given  by 

El  —  aEi 
G\z  =  olGiz 

Gi2  =  OlG\2 

Viz  =  aviz 

Vi2  =  avi2 

Matrix  cracking  is  assumed  if  max{FI(<Ti)}  is  due  to  (T2  or  ere  and  the  corresponding  damage  law 
for  this  mode  is  given  by 

E2  —  OtE2 
G2Z  —  OtG23 
Gi2  =  a<Ji2 
hi  =  OIV21 
V2Z  =  a*'23 

Delamination  is  assumed  if  max{FI{(Ti)}  is  due  to  <73,  (74  or  <75.  The  damage  law  for  this  mode  is 
given  by 

Ez  =  olEz 
G23  —  01G2Z 
Giz  —  olG\z 

hi  =  OLVzi 
Vz2  =  ttl'32 

Appropriate  thermodynamic  constraints  are  applied  to  maintain  a  positive  definite  material  stiff¬ 
ness  matrix  relating  stresses  and  strains. 

Much  research  has  been  devoted  to  the  incorporation  of  microcracking  in  analyses  to  determine 
constitutive  relations.  An  approximate  elasticity  approach  is  presented  in  Reference  [83]  for  de- 
yeloping  reduced  constitutive  relationships  due  to  distributions  of  cracks.  An  idealization  of  a 
composite  laminate  is  made  using  the  unit  cell  depicted  in  Figure  11.  The  unit  cell  is  composed  of 
two  laminae,  lamina  1  is  assumed  to  contain  a  distribution  of  transverse  cracks  while  lamina  2  is 
undamaged. 


(36) 


(35) 


(34) 
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Figure  11.  One-quarter  unit  cell  of  a  composite  laminate. 


The  result  of  the  analysis  yields  reduced  compliance  terms  for  lamina  1  in  laminate  coordinates 
given  by 
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in  which  are  the  undamaged  ply  compliances  in  laminate  coordinates.  The  0i  coefficients  and 


reduced  stiffness  terms,  Q\]\  are  given  by 
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where 
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and  where  /  is  the  crack  spacing  and  are  the  usual  ply  stiffness  coefficients.  An  interesting 
aspect  to  the  developed  expressions  for  reduced  lamina  compliances  is  that  the  compliance  terms 
are  not  only  a  function  of  the  elastic  constants  of  the  cracked  ply  but  are  dependent  on  the  entire 
laminate  construction.  Reference  [84]  presents  a  micromechanically-based  approach  for  determining 
reduced  compliance  terms  due  to  distributed  densities  of  microcracks.  Based  on  the  slip  theory  of 
plasticity,  a  superimposition  of  local  constituent  properties  on  planes  defined  by  individual  crack 
orientations  is  performed  to  determine  macroscopic  properties.  The  resulting  constitutive  model  is 
given  as 

cr  =  Q(¥>)e  (42) 


in  which  (p  is  an  orientation  density  function  of  the  cracking  planes.  The  analysis  is  restricted  to 
2-D  plane  assumptions  of  elasticity  and  to  monotonically  increasing  tensile  loads. 


Laminate  Analysis 


With  the  aim  of  developing  a  numerical  procedure  to  model  general  3-D  composite  laminates,  the 
exact  representation  of  individual  plies  is  considered  impractical.  ‘Exact’  solutions  to  specific  lam¬ 
inate  problems  are  obtained  through  formal  application  of  the  mathematical  theory  of  elasticity 
in  which  powerful  techniques  may  be  applied  to  obtain  solutions  for  a  stratified  or  discontinuous 
elastic  continuum.  These  problems  are  limited  to  highly  idealized  geometries  which  drastically  limit 
their  applicability.  Instead,  the  mechanical  properties  of  an  assembled  laminate  may  be  obtained 
using  a  structural  engineering  approach  known  as  the  Classical  Laminate  Theory  (CLT)  in  which 
an  averaging  of  ply  properties  is  performed  to  obtain  gross  laminate  behavior  [85,  81].  This  ap¬ 
proximation  permits  the  modeling  of  a  larger  class  of  laminate  problems. 
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For  a  3-D  laminate,  lamination  theory  yields  the  material  stiffnesses  coefficients,  C,j,  given  by 


=  («) 

A:=l 

where  H  is  the  total  laminate  thickness,  hk  are  the  individual  ply  thicknesses  and  N  is  the  total 
number  of  plies.  The  transformed  ply  stiffnesses  are  given  by 

On  =  QnCos^O’’  +  Q22^in'^^’‘  +  2Q’l2sin'^e’’cos‘^e'‘  + 

Q22  =  Qnsin‘^d'‘  +  Q22<^os*^'‘  +  2(Q’l2  +  2Q!le)sin‘^9'‘cos‘^e'‘ 

Q44  =  Q’l4CosH'‘  +  Qlsin'^e'‘ 

Qts  =  QLcos^0'‘  +  Qlsin^>‘ 

Q’k  =  Q’k^{Q’{x  +  Q\2-‘^Q\2-^Q’k>i^''9’'<^osH'^  (44) 

On  =  On  +  (Ou  +  O22-2gJ2-4Q^>*ri20Wd'' 

gjg  =  Q\^cos^d’‘  -f- 
^23  =  Q\^sin^6^  -f 

Om  =  (giico^^^^  -  Q’^.iSin^e*^  -  cos2e^{Q\2  +  2Q!^))cos9'‘ sinO^ 

O26  =  -  Q22Cos^9'‘  +  cos29\Q'l2  +  2Ql^))sin9'^ cos9’^ 

O36  =  (0i3  -  Q2z)sin9’^cos9^ 

0$5  =  {Qis-Q\4)sin9'‘cos9'^ 

in  which 

g‘,  =  i;f(i  -vi,4^)HE^Eia’’) 

Q{,  =  Ei(4^  +  4A^)l(E';Eia'‘) 

Qu  =  + ‘'u<'23)/(£i  £fn‘) 

Q22  —  ^2(1  —  ^3^) 

Q33  =  ES(‘'33  +  ‘'M3)/(E!E^a‘‘)  (45) 

=  £5(1  -  4A,)/(EtE^a‘‘) 

Qii  =  g53 

qL  =  G‘, 

ggg  =  g,2 

and 

4,  =  r‘2(£‘/£f) 

‘'L  =  ‘'a(El/El) 

“h  =  <'UEilEi)  (46) 

H*  =  (1  —  *'12^21  —  *'23^32  —  *'l3*'31  “  ^^21^32^13) 

where  Eij,  Gij,  Uij  and  9  are  the  Young’s  moduli,  shear  moduli,  Poisson  ratios  and  fiber  orientation 

angle,  respectively.  The  local  structural  stiffness  is  thus  obtained  using  the  above  averaging  of 
individual  ply  properties.  Specific  problems  may  offer  a  simplification  to  the  above  by  adopting 
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plane-stress  or  plane-strain  as  mptions  in  the  basic  constitutive  relations. 


An  additional  assumption  imi 
however,  plies  tend  to  exhibit 
ing  curing  which  affects  the  Ic 
however,  within  the  order  of  a 
engineering  problems  involvin; 
loading  and  support  conditions 
available  for  use  with  numeri 
lamination  theory  is  a  practic 
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discussion  of  various  lamina  p 
plies  axe  all  /ammafe-depender, 
rections  need  to  be  introduced 
having  to  resort  to  the  nearly 
complex  loading  and  support  ( 


cit  to  CLT  is  that  of  uniform  ply  thickness.  In  real  laminates, 
‘wavy’  distribution  of  thickness  due  to  nonuniform  resin  flow  dur¬ 
al  stress  distribution.  This  local  perturbation  from  uniformity  is, 
aroximation  inherent  to  CLT  and  may  be  neglected.  For  common 
built-up  composite  structures  with  arbitrary  geometries,  complex 
an  approximate  CLT  representation  is  the  only  tractable  approach 
J  techniques  such  as  finite  element  analysis.  Although  classical 
approach  for  simulating  multi-ply  laminates,  the  nature  of  this 
umes  an  inherent  independence  of  individual  plies  which  has  been 
atal  and  analytical  results  to  be  imprecise.  As  noted  above  in  the 
perties,  ply  strength  measures  and  constitutive  laws  for  damaged 
Thus,  the  complexity  of  laminate  analysis  is  ever  present  and  cor- 
t  the  ply  level  to  improve  the  accuracy  of  laminate  theory  without 
tractable  3-D  elasticity  representations  of  general  laminates  under 
nditions. 


Modeling  failure  phenomena  :  the  laminate  level  involves  the  growth  of  local  regions  of  fiber 
breakage  and  matrix  cracking  i  plies  which  coalesce  to  form  significant  3-D  regions  of  damaged 
laminate  material  together  wit  the  extension  of  interply  delaminations  all  of  which  contribute  to 
ultimate  component  failure.  C  mponent  level  failures  are  generally  characterized  as  compressive, 
tensile  or  local/global  buckling  lodes.  The  predominant  failure  modes  are  dependent  on  the  applied 
load  conditions;  under  tensile  ading,  fiber  breakage  and  delamination  growth  causing  significant 
strain  redistribution  are  the  e  ential  modes  involved  in  failure  while  under  compressive  loading 
matrix  cracking,  fiber  breakag  and  local  delamination  buckling  become  the  predominant  failure 
modes.  Simulating  the  gross  b(  avior  of  the  composite  component  in  regards  to  load-displacement 
response  and  local  stress  reco\  ry  may  additionally  involve  nonlinear  material  and  geometric  re¬ 
sponse  under  applied  loads  an(  potential  delamination  growth. 


The  following  subsections  disci  s  laminate-level  behavior  involving  laminate  failure  modes,  nonlin¬ 
ear  load-displacement  response  and  delamination  instability  or  buckling  failure.  It  is  these  effects 
involving  the  gross  laminate  n  ponse  to  applied  loads  and  accumulated  damage  which  must  ulti¬ 
mately  be  simulated  and  asses  d  to  predict  catastrophic  component  failure. 


Laminate  Material  Failui  ' 


A  multitude  of  failure  criteria  for  isotropic  materials  have  been  developed  which  include  single 
component  criteria  such  as  ma  imum  stress,  maximum  strain  and  combined  stress  measures  such 
as  Von  Mises  and  Tresca  based  n  distortional  energy  measures.  None  of  these  theories  differentiate 
between  compression  and  tens  in  and  utilize  a  single  parameter  -  the  yield  strength  -  to  predict 
failure.  The  anisotropic  failurt  criteria  such  as  the  Hoffman,  Tsai-Hill  and  Tsai-Wu  as  discussed 
above  can  be  applied  in  an  app  oximate  sense  utilizing  transformed  and  averaged  ply  strengths  to¬ 
gether  with  laminate  stresses,  'onservative  approximations  of  local  laminate  failure  have  utilized 
the  first  ply  failure  (FPF)  critc  ia  for  predicting  local/global  laminate  failure.  These  criteria,  how- 
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ever,  are  applied  using  the  state  of  stress /strain  at  a  point.  For  the  prediction  of  complete  laminate 
failure,  global  criteria  are  required  to  assess  the  integrity  of  composite  components.  Local  regions 
of  material  damage  cause  a  general  redirection  of  load  ‘flow’  around  regions  of  reduced  stiffness  as 
depicted  in  Figure  12. 


Figure  12.  Load  redistribution  due  to  local  material  stiffness  loss. 


The  redistribution  of  load  thereby  intesifies  the  state  of  stress  in  adjacent  areas  and  tends  to 
precipitate  additional  material  failure  or  enhance  delamination  growth  leading  to  buckling  failure. 

Laminate  level  failure  criteria  may,  therefore,  be  based  on  gross  measures  of  maximum  strains 
under  applied  tensile  or  compressive  loading  or  reduction  in  flexural  stiffness.  In  an  analytical  sim¬ 
ulation,  significant  reductions  in  stiffness  along  any  cross-sectional  plane  could  be  used  to  indicate 
catastrophic  failure.  Alternatively,  distributed  damage  resulting  in  excessive  strain  or  a  drop  in 
the  integrated  resistance  to  applied  loading  throughout  any  cross  section  may  be  considered  as 
indicative  of  ultimate  component  failure. 


Nonlinear  Laminate  Behavior 

Nonlinear  laminate  behavior  may  be  manifest  in  geometric  nonlinearity  involving  large  displace¬ 
ments  and  buckling  or  in  material  behavior  involving  nonlinear  elastic,  viscoelastic  and/or  plastic 
response  to  applied  loads.  The  current  study  neglects  inelastic  ply  material  response  and  assumes, 
at  most,  a  nonlinearily  elastic  behavior  to  failure.  Ply  level  nonlinear  material  behavior  dictate 
the  global  elastic  behavior  of  the  composite  component.  As  mentioned  in  Section  3.2,  nonlinearity 
is  primarily  observed  in  the  ply  shear  modulus.  In  averaging  individual  ply  properties,  general 
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cross-ply  laminates  such  as  [±45]  layups  demonstrate  significant  nonlinear  laminate  response  while 
quasi-isotropic  layups  such  as  [0/  ±  45/90]  behave  linearily.  Thus,  some  discretion  is  available  to 
the  analyst  in  the  choice  of  linear  or  nonlinear  representation  of  gross  laminate  behavior  based  on 
laminate  ply  composition.  For  arbitrary  laminates,  stress  analyses  must  account  for  the  nonlinear 
load-displacement  behavior  of  the  laminate  which  requires  an  iterative  solution  scheme  for  all  but 
the  simplest  geometries.  In  addition,  large  displacements  constitute  a  nonlinear  geometric  response 
which  must  be  accounted  for  in  analytical  simulation  by  including  nonlinear  terms  in  the  stress- 
strain  relations.  In  real  materials,  both  nonlinear  material  and  geometric  responses  are  generally 
significant  in  predicting  the  response  to  extreme  applied  loading  in  order  to  accurately  calculate 
stresses  for  residual/failure  strength  prediction.  Analytic  or  numerical  modeling  approaches  are, 
therefore,  tremendously  complicated  by  these  nonlinear  effects  in  obtaining  accurate  solutions. 


Laminate  Delamination:  Initiation,  Buckling  and  Stability 


Under  compressive  loading,  debonding  of  adjacent  plies  or  delaminations  may  be  considered  as 
the  most  severe  form  of  laminate  damage  in  contributing  to  the  reduction  in  residual  strength 
through  local  buckling  failure  and  flexural  stiffness  reduction.  In  addition  to  the  presence  of  initial 
delaminations  due  to  manufacturing  defects  or  impact  events,  laminate  stresses  can  initiate  and 
propagate  new  delaminations.  Numerous  criteria  have  been  developed  to  predict  the  initiation  of 
delaminations  [86,  87,  88].  For  example.  Reference  [89]  presents  a  quadratic  delamination  criterion 
of  the  following  form 


where  all  stresses  are  average  measures  computed  as 

1  fXave 

aij  =  — —  /  cTijdx  (48) 

-A  at/e  Jo 

where  Xave  is  a  characteristic  length  either  assumed  or  determined  experimentally.  The  strength 
parameters  are  given  as 

Zt  =  tensile  interlaminar  normal  strength. 

Zc  =  compressive  interlaminar  normal  strength. 

Zai  =  interlaminar  shear  strength  for  Cis  stresses. 

Zs2  =  interlaminar  shear  strength  for  <T23  stresses. 

Local  buckling  of  laminate  sublayers,  as  depicted  in  Figure  13,  can  create  high  interlaminar  stresses 
at  the  delamination  fronts  which  may  cause  further  extention  of  the  buckled  layer.  Consequentially, 
the  redistribution  of  loads  may  induce  additional  delaminations  and  inplane  material  damage. 
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Figure  13.  Local  buckling  of  sublaminate  due  to  delamination. 


Delaminations  may  be  formed  during  laminate  construction,  through  microcrack  coalescence  un¬ 
der  service  fatigue  loading  or  through  impact  with  foreign  objects.  In  built-up  laminates  it  has 
been  noted  that  delamination  due  to  impact  tends  not  to  occur  between  plies  of  like  orientation 
[90,  8].  Various  analytical  models  have  been  developed  using  structural  engineering  approximations 
such  as  beams  and  plates  to  represent  laminate  sublayers  [6,  10,  7,  11],  and  approximate  elasticity 
solution  approaches  have  been  formulated  in  [9,  11]  using  energy  methods.  Other  investigations 
have  utilized  specialized  finite  element  analyses  which  have  determined  that  the  stress  state  around 
delaminated  regions  is  generally  biaxial  and  that,  depending  on  ply  layup  and  delamination  con¬ 
figuration,  buckling  can  initiate  under  applied  tensile  loading  [6,  11,  12].  Additional  eflFects  which 
have  been  studied  involve  bridging  between  delamination  surfaces  due  to  fibers  or  stitching  repair 
[91]. 


Buckling  or  structural  instability  may  be  defined  as  the  condition  at  which  the  change  in  the  total 
potential  energy  associated  with  an  arbitrary  displacement,  6q,  from  equilibrium  is  an  extremum 
or,  equivalently,  that  the  total  potential  satisfies  the  Trefftz  criterion  [92]  given  by 


n  =  u-i-v 


(49) 


and 


(50) 


in  which  U  is  the  total  strain  energy,  V  is  the  potential  energy  and  (f>  is  the  generalized  displacement. 
Crack  stability  leading  to  delamination  growth  has  been  most  successfuly  investigated  using  a 
fracture  mechanics  approach  [11,  93,  94].  The  basic  Griffith  energy  criteria  for  crack  growth  may 
be  stated  as 

|(W-U)  =  G.  (51) 

where  W  is  the  work  performed  by  external  applied  loading,  U  is  the  internal  elastic  strain  energy, 
Gc  is  the  critical  strain  energy  release  rate  required  to  create  new  fracture  surface  and  a  is  the 
initial  crack  length.  A  versatile  technique  for  calculating  the  strain  energy  release  rate  with  ideal 
finite  element  suitability  is  the  crack-closure  integral  method  [95,  96].  The  essential  idea  is  that  the 


energy  absorbed  in  extending  a  crack  by  an  amount  Aa  is  equal  to  the  work  required  to  close  the 
crack  to  its  ori^nal  length.  The  energy  statement  may  be  expressed  as 

G  =  Gi  +  Gii  (52) 

where  G  is  the  total  strain  energy  release  rate  and  G/  and  G//  are  the  energy  release  rates 
associated  with  Mode  I  and  Mode  II  deformation  components.  With  the  polar  coordinate  system 
depicted  in  Figure  14,  the  strain  energy  release  rate  components  may  be  expressed  as 


G/  =  lim  — ^ 
Aa-*0  2Ao , 

/  (Ty(Aa  -  r,  0)i)(r,  7r)dr 

Jo 

(53) 

G//  =  lim 

Aa— 0  2Aa , 

rAa 

/  TxyiAa  -  r,  0)u(r,  Tr)dr 

Jo 

(54) 

where  Oy  and  r^y  are  stresses  near  the  crack  tip  -  or  delamination  front,  u  and  v  are  the  relative 
Mode  I  sliding  and  Mode  II  opening  displacements  between  opposing  points  along  the  crack  length 
and  Aa  is  the  incremental  extension  in  the  initial  crack  length,  a. 


Figure  14.  Coordinate  system  at  crack  front 


The  presence  of  an  initial  profile  of  delaminations  through  the  laminate  thickness  may,  thus,  exhibit 
growth  during  subsequently  applied  loads  due  to  coupling  between  buckling  and  crack  instability 
leading  to  component  failure.  In  purely  structural  instability  situations  in  which  growth  of  initial 
delaminations  is  not  predicted,  local  buckling  may  couple  with  global  buckling  to  cause  failure 
through  a  ‘mixed-mode’  instability. 


Finite  Element  Modelling 

Numerous  specialized  analytical  approaches  to  predicting  residual  strength  in  composite  materials 
have  been  made  incorporating  specific  failure  modes  and  highly  idealized  geometric  configurations. 
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The  finite  element  method  provides  a  completely  generic  numerical  approach  to  potentially  incor¬ 
porate  all  salient  material  behavior  and  failure  modes  while  allowing  arbitrary  geometries,  applied 
loads  and  support  conditions  to  be  incorporated  in  the  analysis  [97].  The  present  study  is  ultimately 
directed  towards  simulating  impact  damage  in  thick  laminated  composites  which  immediately  sug¬ 
gest  plate-type  finite  element  formulations  to  represent  the  elastic  continuum.  In  the  domain  of 
single-layer  representations  simple  first-order  shear  deformable  theories  such  as  presented  in  Ref¬ 
erence  [98]  or  higher-order  formulations  such  as  those  presented  in  References  [100-103]  may  be 
considered.  However,  an  important  consideration  is  the  representation  of  three  dimensional  effects 
in  any  plate  formulation.  Reference  [103]  presents  a  viable  equivalent  single-layer  plate  theory  incor¬ 
porating  three-dimensional  behavior  while  requiring  exclusively  Poisson-type  boundary  conditions. 
Layer-wise  theories  have  been  developed  such  as  those  of  References  [104,  105]  which  attempt  to 
account  in  greater  detail  for  variations  through  the  laminate  thickness.  One  drawback,  however, 
inherent  to  all  plate  formulations  is  the  impossibility  of  representing  delamination  growth  and  local 
sublaminate  buckling  which  may  be  a  dominant  failure  mode  in  certain  laminate  types  and  damage 
distributions.  Therefore,  for  maximum  versatility  in  three  dimensional  modeling  of  thick  laminates, 
an  8-node  solid  brick  element  is  suggested  which  utilizes  the  hybrid  stress  technique  to  represent 
elastic  stiffness  response  while  purely  displacement-based  field  assumptions  are  used  to  develop 
differential  stiffness  relations  to  model  geometric  nonlinear  response.  Material  degradation  due  to 
sequential  failures  is  determined  using  the  above  failure  criteria  together  with  assumed  ply  damage 
laws  and  is  directly  incorporated  into  the  element  formulation.  The  basic  element  formulations 
together  with  nonlinear  material  and  geometric  solution  algorithms  are  presented  in  the  following 
subsections. 


Elastic  Stiffness  Formulation 

The  hybrid  stress  technique  is  utilized  to  form  elastic  stiffness  relationships.  Details  on  the  hybrid 
method  may  be  found  in  References  [108-114].  The  structure  of  element  matrices  are  defined  by 
the  energy  functional  which  is  selected  as  the  extended  Hellinger-Reissner  functional  and  may  be 
stated  as 


Hr  =  j[{-l/2)a^Sfr  +  a^{Luq)-{L^a)^ux]dv  (55) 

where  a  is  the  assumed  stress  field,  S  is  the  material  compliance  matrix,  u,  and  ua  are  the  assumed 
compatible  and  incompatible  displacement  fields,  and  L  is  the  differential  operator  relating  strains 
to  displacements. 

The  assumed  stresses  may  be  represented  by 


a  =  P(3  (56) 

where  P  is  a  matrix  of  polynomial  terms  and  is  a  vector  of  undetermined  expansion  coefficients. 
The  displacement  field  is  assumed  over  the  element  domain  as 

u  =  u,  +  UA  (57) 

in  which  compatible  and  incompatible  displacement  components  are  given  by 

u,  =  Nq 

ua  =  MA 
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(58) 

(59) 


where  N  and  M  are  compatible  and  incompatible  displacement  shape  functions,  respectively,  q 
axe  nodal  displacements,  and  A  are  Lagrange  multipliers  which  enforce  internal  constraints.  In  the 
form  of  (55),  performing  the  variation  with  respect  to  ua,  the  incompatible  displacements  may  be 
used  to  vaxiationally  enforce  a  priori  the  field  equilibrium  conditions  through  the  last  term  in  (55) 
which  requires  that 

S  J  (L^  a)^  uxdv  =  0  (60) 

or 

L^a  =  0  (61) 

Applying  the  constraints  to  the  stress  modes  results  in  the  reduced  functional 

Ur=  f  [(-l/2)<T^Scr  +  o-^(Lu5)]du  (62) 


Substituting  (56),  (58)  and  (59)  into  (55)  yields 


Ur  =  J  [(-l/2)/3^P^SP/3  +  /3^P^(LN)q]du 

(63) 

or 

n«  =  (-l/2)/3^H/3  +  ;3^Gq 

(64) 

where 

H  =  iP^SPdv 

(65) 

G  =  /„P^(LN)dn  =  /„P^Bdn 

(66) 

Seeking  a  stationary  value  of  the  functional  by  taking  the  first  variation  with  respect  to  /3  yields 


(3  =  H-^Gq  (67) 

substituting  the  resulting  expression  for  (3  into  (64),  the  variation  with  respect  to  q  yields  the 
element  stiffness  matrix  as 

K  =  G^H-^G  (68) 

An  8-node  solid  element  as  depicted  in  Figure  15  is  selected  to  represent  the  laminate  continuum. 


Figure  15.  Hexahedral  element  configuration. 
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The  displacement  functions  u,  are  ^ven  by 


u 


9 


1 

*  1 

r, 

^  =  E  o  (1  +  ^.0(1  +  ’?.»7)(i  +  C.C)  ^ 

t=l  ® 

Ui  1 
Vi  \ 


(69) 


The  isoparametric  mapping  between  physical  and  natural  coordinates  is  given  by 
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The  stress  field  is  initially  assumed  as  complete  quadratic  expansions  in  natural  or  tensor  coordi 
nates  given  by 


where 


[Po] 


[Po] 


p  = 


[Po]  = 


(72) 


(73) 


As  detailed  in  Reference  [110],  incompatible  displacement  functions  are  assumed  such  that  the 
polynomial  order  the  total  displacement  field  is  cubic  thereby  yielding  a  resultant  quadratic  strain 
field  of  the  same  order  as  the  assumed  stresses.  The  assumed  stress  modes  are  then  subjected  to 
the  constraint  given  by  equation  (61)  yielding  an  element  stress  field  with  18  independent  stress 
modes  given  by 
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(74) 


To  preserve  the  constant  stress  terms,  the  natural  stresses  are  mapped  to  physical  space  through  a 
contravariant  transformation  using  Jacobians  computed  at  the  element  centroid 


=  WHJoYjr 

The  stress  field  expressed  in  physical  or  Cartesian  coordinates  is  given  by 
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The  above  8-node  solid  element  formulation  based  on  hybrid  element  theory  is  a  robust  element 
with  excellent  behavior  in  normal  and  bending  deformation  modes. 


Geometric  Stiffness  Formulation 


With  the  discretized  geometry  inherent  to  the  finite  element  technique,  the  effects  of  nonlinear 
strains  may  be  represented  by  discrete  matrix  formulations  based  on  the  current  stress  state.  The 
geometric  stiffness  matrix  is  formulated  using  the  assumed  displacement  fields  in  (69)  and  is  given  by 

[K,]  =  J[BNLf[T][BNL]dv  (80) 


where 
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and 


The  Bnl  matrix  is  given  by 


where 


(^xxO 

“^xyO 

[S]  = 

TxyO 

<^yyO 

TyzO 

'^yxO 

OzzQ 

[Ba^l]  = 

■[G] 

0 

0 

[G] 

0 

0 

0 

0 

[G] 

Ni^ 

0 

0 

N2^ 

...  Ng^ 

[G]  = 

Ni,v 

0 

0 

N2.v 

-  N8,„ 

Ni,. 

0 

0 

N2,, 

...  Ng,, 

(82) 


(83) 


(84) 


The  modeling  of  nonlinear  strain  response  may  be  extended  to  include  large  strain  and  displacement 
behavior  through  selection  of  standard  algorithmic  solution  procedures  in  the  finite  element  method 
[113,  115].  These  methods  involve  iterative  procedures  utilizing  total  or  updated  Lagrangian  rep¬ 
resentations  of  the  structural  model  kinematics. 
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Delamination  Propagation 


As  discussed  above  in  section  4.3,  the  crack-closure  integral  method  [95,  96]  may  be  used  for  cal¬ 
culating  the  strain  energy  associated  with  potential  delamination  propagation.  A  particularily 
attractive  feature  of  this  technique  is  the  ability  to  use  a  coarse  mesh  without  recourse  to  singu¬ 
larity  elements  to  obtain  accurate  measures  of  the  strain  energy  release  rate.  Within  a  standard 
finite  element  approximation  framework,  the  evaluation  of  strain  energy  release  rates  is  based  on 
nodal  forces  and  displacements.  The  delamination  is  represented  by  a  surface  of  coincident  nodes 
which  are  used  to  create  a  plane  of  discontinuity  through  the  definition  of  element  connectivities. 
A  representative  finite  element  discretization  in  the  region  of  the  crack  tip  is  shown  in  Fiqure  16. 
The  work  required  to  close  the  incremental  growth  in  the  crack  length  may  be  expressed  as 
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Figure  16.  Finite  element  model  of  delamination  front. 


In  the  computation  of  the  x-  and  y-  nodal  forces  required  to  hold  the  nodes  at  c  and  d  together,  the 
values  for  Fc  and  Tc  are  approximated  by  the  nodal  forces  at  nodes  e  and  /  assuming  an  inverse 
square  root  singular  behavior  near  the  crack  tip  as 

F.  =  (8T) 

*2 

and 

Tc  =  (88) 

*2 

Delamination  growth  is  thus  predicted  when  the  release  of  elastic  energy  associated  with  the  ex¬ 
tension  of  the  delamination  exceeds  the  critical  strain  energy  release  rate  of  the  material  such 
that 

G  =  Gi  +  Gn>Gc  (89) 

When  the  above  criterion  is  met,  the  local  delamination  front  is  extended  at  the  evaluation  node 
and  incorporated  into  the  model. 

The  accurate  prediction  of  delamination  growth  requires  a  local  refinement  of  the  finite  element 
model  at  the  delamination  front.  Predicting  sequential  growth  at  potentially  multiple  delamina¬ 
tions  implies  a  high  degree  of  mesh  refinement  and  iteration  in  which  the  connectivity  of  the  finite 
element  model  is  updated  to  account  for  delamination  extension.  The  crack  closure  technique  does 
permit  a  degree  of  coarseness  at  the  crack  tip  as  demonstrated  in  Reference  [95]  in  which  ratios  of 
incremental  crack  extension,  Ac,  to  total  crack  length,  c,  up  to  20%  yielded  inaccuracies  less  than 
6%.  Further  study  is  required  to  extend  the  above  crack  closure  methodology  to  bimaterial  inter¬ 
face  problems  which  possess  a  general  A  -  order  singularity.  Additional  investigation  is  warrented 
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to  compare  the  computational  requirements  with  accuracy  in  predicting  sequential  delamination 
growth  in  multiply  delaminated  laminates. 


Material  Nonlinearity  Solution  Algorithm 


Assuming  infinitesimal  displacements  and  small  strains,  at  some  applied  load  level,  Pa,  the  stress/strain 
relationship  may  be  expressed  as 

=  C*(e)t*  (90) 

where  C’‘{e)  is  the  strain  dependent  material  stiffness  matrix,  cr*  and  are  the  stress  and  strain 
vectors  at  the  load  level.  To  solve  the  nonlinear  global  nonlinear  system  of  equations,  a  Newton- 
Raphson  sheme  may  be  employed  in  which  an  iterative  approach  is  used  to  obtain  the  equilibrium 
state 

K.  AD,-+i  =  Pa  -  Ri  (91) 

Di+i  =  Di  +  ADi+i  (92) 

where  K,-  is  the  elastic  stiffness  matrix  at  the  iteration.  Pa  are  the  applied  loads,  are  the 
current  internal  recovered  loads  and  AD,+i  are  the  incremental  displacements  which  are  summed 
to  give  the  total  displacements  for  the  iteration.  This  scheme  is  depicted  in  Figure  17. 
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Figure  17.  Newton-Raphson  iteration  scheme. 


Iterations  are  continued  until  the  difference  between  the  vector  norms  of  the  externally  applied 
loads  and  the  internal  recovered  loads  are  less  than  a  prescribed  tolerance 

|P„  -  R.|  <  Tol  (93) 

Convergence  rates  may  be  enhanced  by  performing  line  searches  in  the  current  ‘Newton’  direction 
defined  by  AD,.|.i  such  that 

Dj+i  =  Dj  +  /3AD,+i  (94) 

where  (3  is  determined  through  line  search  methods  such  as  bisection,  secant,  Ridder  or  Regula  Falsi 
[113,  115,  116]. 


Geometric  Nonlinear  Solution  Algorithm 

In  the  linear  theory  of  buckling,  a  solution  is  sought  at  the  point  where  the  locus  of  equilibrium 
states  bifurcate  as  depicted  in  Figure  18. 


DISPLACEMENT 

Figure  18.  Bifurcation  of  equilibrium  states. 

Buckling  is,  therefore,  a  nonlinear  geometric  phenomemon.  In  order  to  approximate  this  behavior, 
higher-order  strain  terms  -  usually  quadratic  -  are  used  to  compute  the  differential  stiffness  matrix 
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K(7  as  discussed  in  Section  5.2  resulting  in  the  nonlinear  equilibrium  relationship 


1 


(K  +  K^)X=R  (95) 

To  determine  the  load  level  at  which  bifurcation  occurs  requires  the  solution  to  the  following  general 
eigenvalue  problem 

(K  +  iK^)X  =  0  (96) 

where  the  critical  buckling  load  is  then  given  by 

Rcrit  =  AR  (97) 

In  order  to  determine  the  lowest  eigenvalues  of  the  global  system,  a  subspace  iteration  is  employed 
[113,  114].  An  initial  set  of  iteration  vectors,  Xfc,  are  selected  to  span  the  space  defined  by  the 
desired  eigenvectors,  Iterations  are  then  performed  to  converge  the  eigenvectors  and  eigenvalues 
in  the  limit 

X/t  — >  ^  as  k  — >•  oo  (98) 

The  subspace  mapping  is  defined  by 

KXfc+i  =  K,Xfc  (99) 

The  mapping  or  iteration  vectors  are  selected  in  accordance  with  a  procedure  outlined  in  Reference 
[113].  The  selection  is  specified  such  that  the  first  column  in  X  is  taken  as  the  diagonal  in  Ka  while 
remaining  columns  are  assigned  a  unit  value  in  the  row  position  corresponding  to  min{kiifkaii)- 
The  global  elastic  and  differential  stiffness  matrices  are  then  mapped  into  the  subspace  as 

Kfc+i  =  Xf+iKXfc+i  (100) 

K<7fc+1  =  (101) 

resulting  in  the  reduced  eigenvalue  problem  given  by 

(102) 

where  A  is  a  diagonal  matrix  of  the  system  eigenvalues.  Since  the  initial  iteration  vectors  are 
usually  approximate,  an  improvement  to  the  system  eigenvectors  is  given  by 

Xjfe+i  =  Xk+iQk+i  (103) 


and  an  iteration  of  the  subspace  is  continued  until  the  eigenvalues  reach  a  stationary  value  defined 
by  a  tolerance  as 


-  A.*| 
|Af+'| 


<tol  ,  i  =  1,2,3,  ...,n 


(104) 


During  iterations  a  check  may  be  made  to  ensure  that  the  smallest  eigenvalues  are  being  obtained 
by  using  the  Sturm  sequence  property  of  the  characteristic  polynomials  of  eigensystems.  If  a  shift, 
fi,  is  introduced  slightly  greater  than  the  current  estimate  of  A„,  the  factorization 


K,  =  K  -  =  LDL^ 


(105) 


yields  a  diagonal  matrix,  D,  in  which  the  number  of  negative  diagonals  is  equal  to  the  number  of 
eigenvalues  smaller  than  /i. 
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An  important  consideration  in  modeling  impact  damaged  composites  is  the  presence  of  multi¬ 
ple  delaminations  schematically  depicted  in  Figure  19. 


Figure  19.  Laminate  with  multiple  embedded  delaminations. 

The  stability  analysis  of  multiply  delaminated  laminates  is  complicated  by  the  existence  of  a  number 
of  potential  instability  modes.  If  delaminations  are  simply  modeled  as  planes  of  element  discon¬ 
tinuity  and  a  linear  buckling  analysis  as  detailed  above  is  applied,  predicted  buckling  modes  will 
include  physical  modes  and  nonphysical  modes  in  which  the  displacement  of  interior  sublayers  de¬ 
form  through  surrounding  laminate  material.  This  behavior  is  depicted  in  Figure  20. 


P  '=!> 


<1^P 
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The  investigation  of  buckling  in  laminates  containing  multiple  delaminations  has  not  been  ade- 
quetely  addressed  in  the  literature.  Reference  [12]  considers  a  finite  element  study  of  symmetrically 
located  delaminations  and  uses  the  crack-closure  technique  to  suppress  nonphysical  or  inadmiss- 
able  modes.  It  is  noted  therein  that  when  one  delamination  is  undergoing  buckling,  the  presence 
of  an  additional  delamination  represented  as  a  plane  of  frictionless  contact  has  only  a  slight  effect 
»  on  predicted  buckling  loads  over  those  predicted  with  the  second  delamination  removed  from  the 

model.  The  buckling  and  growth  of  multiple  delaminations  is  studied  in  References  [117,  118] 
using  Von  Karman  plate  theory  in  a  finite  element  context  to  model  the  sublayers.  An  iterative 
''  scheme  using  the  conjugate  gradient  method  is  utilized  to  determine  equivalent  contact  pressures 

to  avoid  nonphysical  buckling  modes.  Additional  work  on  the  imposition  of  contact  constraints  in 
finite  element  buckling  analysis  may  be  found  in  References  [119,  120].  The  analytical  approaches 
demonstrate  the  importance  of  sublayer  contact  in  the  calculation  of  strain  energy  release  rates 
and  may  provide  a  basis  for  further  development  to  accurately  model  multiple  delaminations  in 
composites. 


Combined  Analysis  for  Residual  Strength  Prediction 


The  combined  effects  of  accumulated  material  damage,  inherent  material  nonlinearity,  and  po¬ 
tential  delamination  growth  and  sublaminate  buckling  significantly  complicates  a  residual  strength 
analysis.  The  material  nonlinearity  together  with  a  stability  analysis  requires  a  piece-wise  linear 
approach  in  which  applied  loads  are  incremented  and  eigenvalues  extracted  with  buckling  predicted 
when  the  largest  eigenvalue  is  unity.  In  performing  such  a  procedure,  each  iteration  requires  subit¬ 
erations  to  converge  the  prediction  of  ply  failures  and  delamination  growth  because  each  alters  the 
local  stiffness  distribution  in  the  laminate  thereby  upsetting  equilibrium  and  internal  load  distri¬ 
bution.  A  schematic  of  the  complete  solution  procedure  is  depicted  in  Figure  21. 


39 


INPUT  INITIAL  DAMAGE  STATE  AND  BEGIN 
NEWTON-RAPHSON  ITERATION 


[K]n{ADi^i}  =  {ARi+i} 
{Di+i}  =  {Di}  +  {ADi} 


TEST  FOR  ADDITIONAL  DAMAGE 

MATRIX  CRACKING:  >  1? 
FIBER  BREAKAGE;  ef>l?  _ 
CRACK  GROWTH:  Aa  >  0? 
BUCKLING  FAILURE:  Xmax=1? 

<NO> 

T _ 

-<NO>— ]  |{Pn^.iWRi..i}|<TOL  - 

<YES> 

r 

- {Pn+i}  =  {Pn}  +  {AP} 


<NO> 


INCORPORATE 

DAMAGE 


COMPLETE 

LAMINATE 

FAILURE? 

<YES> 

i _ 

OUTPUT 

RESIDUAL 

STRENGTH 


Figure  21,  Flow  chart  of  combined  iterative  solution  algorithm. 


Closure 


Various  modeling  considerations  in  the  prediction  of  residual  strength  in  built-up  composite  lami¬ 
nates  have  been  discussed.  Regardless  of  the  inherent  accuracy  of  analytic  modeling,  the  prediction 
of  residual  strength  is  strongly  dependent  on  the  accuracy  of  the  inputed  initial  damage  state.  The 
three-dimensional  state  of  damage  may  be  obtained  through  various  NDE  technologies,  such  as 
Acoustic/Ultrasonics,  Thermography  and  Radiography,  each  yielding  different  degrees  of  resolu¬ 
tion  in  the  detection  of  various  internal  damage  modes. 


The  mathematical  representation  of  composites  must  deal  with  three  scale  regimes:  The  micro¬ 
scopic,  the  ply  and  the  composite  component.  In  each  regime,  approximations  are  introduced  to 
make  the  analysis  tractable.  At  the  microscopic  level,  the  characteristics  of  fiber  and  matrix  fail¬ 
ure  may  be  analyzed  through  micromechanical  models  to  develop  failure  criteria  and  to  deduce 
macroscopic  behavior.  At  the  ply  level,  the  assumptions  include  the  representation  of  a  lamina  as 
a  homogeneous  elastic  continuum,  the  adoption  of  ‘independent’  or  ‘interactive’  ply  failure  modes, 
degradation  of  elastic  moduli  through  developed  damage  laws,  and  the  approximation  of  nonlin¬ 
ear  lamina  behavior.  At  the  component  level,  approximations  are  made  regarding  the  averaging 
of  ply  material  properties  to  develop  continuous  laminate  elastic  constants,  delamination  growth 
and  sublaminate  buckling  behavior.  Additional  approximations  are  introduced  through  simplifying 
assumptions  of  component  geometry,  applied  loading  and  support  conditions.  The  variablity  in  the 
construction,  properties  and  defects  in  composite  components  contribute  to  the  usually  wide  scat¬ 
ter  observed  in  experimental  tests  from  specimen  to  specimen  and  present  an  inherent  statistical 
aspect  to  any  attempt  to  mathematically  simulate  actual  composite  behavior  to  predict  residual 
strength. 

With  a  consideration  of  the  aforementioned  caveats,  a  finite  element-based  numerical  approach 
has  been  detailed  to  determine  the  maximum  load  carrying  capability  in  damaged  laminates  by 
outlining  how  various  laminate  phenomena  can  be  incorporated.  The  versatility  of  the  finite  ele¬ 
ment  method  in  representing  arbitrary  geometries,  boundary  conditions  and  load  conditions  while 
accounting  for  material  failure  together  with  nonlinear  behavior  and  delamination  growth  leading 
to  local  buckling  instability  should  provide  a  more  accurate  representation  of  all  relevant  failure 
mechanisms  exhibited  in  real  laminates.  These  phenomena  are  incorporated  through  solution  algo¬ 
rithms  such  as  Newton-Raphson  for  material  nonlinearity  and  subspace  iteration  for  the  prediction 
of  buckling  failure.  Additional  considerations  such  as  the  incorporation  of  large  displacement  and 
large  strain  effects  in  simulating  the  behavior  of  laminates  after  initial  failure  can  be  made  through 
standard  procedures  existent  in  finite  element  theory. 


■> 


41 


References 

[1]  Madaras,  E.  I.,  Poe,  C.  C.  and  Heyman,  J.  S.,  “Combining  Fracture  Mechanics  and  Ultra¬ 
sonic  NDE  to  Predict  the  Strength  Remaining  in  Thick  Composites  Subjected  to  Low-Level 
Impact,”  1986  Ultrasonics  Symposium,  pp.  1051-1059. 

[2]  Garrett,  R.  A.,  “Effect  of  Manufacturing  Defects  and  Service-Induced  Damage  on  the  Strength 
of  Aircraft  Composite  Structures,”  Composite  Materials:  Testing  and  Design  (Seventh  Con¬ 
ference),  ASTM  STP  893,  American  Society  for  Testing  and  Materials,  pp.  5-33,  (1986). 

[3]  Cairns,  D.  S.  and  Lagace,  P.  A.,  “Residual  Tensile  Strength  of  Graphite/Epoxy  and  « 

Kevlar/Epoxy  Laminates  with  Impact  Damage,”  Composite  Materials:  Testing  and  Design 

(Ninth  Volume),  ASTM  STP  1059,  American  Society  for  Testing  and  Materials,  pp.  48-63, 

(1990). 

[4]  Caprino,  G.,  “Residual  Strength  Predictions  of  Impacted  CFRP  Laminates,”  J.  Comp.  Mats., 

18,  pp.  508-518,  (1984). 

[5]  Lai,  K.  M.,  “Prediction  of  Residual  Tensile  Strength  of  Transversely  Impacted  Composite 
Laiiiinates,”  Research  in  Structural  and  Solid  Mechanics  -  1982,  NASA  Conference  PubUcar 
tion  2245,  pp.  97-112. 

[6]  Shu,  D.  and  Mai,  Y.  W.,  “Buckling  of  delaminated  composites  re-examined,”  Comp.  Sci. 

Tech.,  47,  pp.  35-41,  (1993). 

[7]  Chai,  H.,  Babcock,  C.  D.  and  Knauss,  W.  G.,  “One  Dimensional  Modelling  of  Failure  in 
Laminated  Plates  by  Delamination  Buckling,”  Int.  J.  Solids  Struct.,  Vol.  17,  No.  11,  pp.l069- 
1083, (1981). 

[8]  Lui,  D.,  ‘Delamination  Resistance  in  Stitched  and  Unstitched  Composite  Plates  Subjected  to 
Impact  Loading,’  J.  Reinf.  Plastics  Compos.,  9,  pp.  59-69,  (1990). 

[9]  Simitses,  G.  J.,  Sallam,  S.  and  Yin,  W.  L.,  “Effect  of  Delamination  of  Axially  Loaded  Homo¬ 
geneous  Laminated  Plates,”  AIAA,  Vol.  23,  No.  9,  pp.  1437-1444,  (1985). 

[10]  Shivakumar,  K.  N.  and  Whitcomb,  J.  D.,  “Buckling  of  a  Sublaminate  in  a  Quasi-Isotropic 
Composite  Laminate,”  J.  Comp.  Mats.,  19,  pp.  2-18,  (1985). 

[11]  Wang,  S.  S.,  Zahlan,  N.  M.  and  Suemasu,  H.,  “Compressive  Stability  of  Delaminated  Random 

Short-Fiber  Composites,  Part  I  -  Modeling  and  Methods  of  Analysis,”  J.  Comp.  Mats.,  19,  < 

pp.  296-316,  (1985). 

[12]  Wang,  S.  S.,  Zahlan,  N.  M.  and  Suemasu,  H.,  “Compressive  Stability  of  Delaminated  Random  ^ 

Short-Fiber  Composites,  Part  II  -  Experimental  and  Analytical  Results,”  J.  Comp.  Mats., 

19,  pp.  317-333,  (1985). 

[13]  Madenci,  E.  and  Westmann,  R.  A.,  “Local  Delamination  Growth  in  Layered  Systems  Under 
Compressive  Load,”  J.  Appl.  Mech.  ASME,  60,  pp.  895-902,  (1993). 

[14]  Madenci,  E.  and  Westmann,  R.  A.,  “Local  Delamination  Buckling  in  Layered  Systems,”  J. 

Appl.  Mech.  ASME,  58,  pp.  157-166,  (1991). 


42 


[15]  Chang,  F.  K.  and  Chang,  K.  Y.,  “Post- Failure  Analysis  of  Bolted  Composite  Joints  in  Tension 
or  Shear-Out  Mode  Failure,”  J.  Comp.  Mats.,  Vol.  21,  pp.  809-833,  (1987). 

[16]  Chang,  F.  K.,  Scott,  R.  A.  and  Springer,  G.  S.,  “Failure  Strength  of  Nonlinearly  Elastic 
Composite  Laminates  Containing  a  Pin  Loaded  Hole,”  J.  Comp.  Mats.,  18,  pp.  464-477, 
(1984). 

[17]  Chang,  F.  K.  and  Chang,  K.  Y.,  “A  Progressive  Damage  Model  for  Laminated  Composites 
Containing  Stress  Concentrations,”  J.  Comp.  Mats.,  Vol.  21,  pp.  834-855,  (1987). 

[18]  Tian,  Z.  and  Swanson,  S.  R.,  “Residual  Tensile  Strength  Prediction  on  a  Ply-by-Ply  Basis 
for  Laminates  Containing  Impact  Damage,”  J.  Comp.  Mats.,  26,  pp.  1193-1206,  (1992). 

[19]  Reddy,  Y.  S.  and  Reddy,  J.  N.,  “Three-Dimensional  Finite  Element  Progressive  Failure  Analy¬ 
sis  of  Composite  Laminates  Under  Axial  Extension,”  J.  of  Composites  Technology  &  Research, 
Vol  15,  No.  2,  pp.  73-87,  (1993). 

[20]  Reddy,  J.  N.,  “On  Refined  Computational  Models  of  Composite  Materials,”  Int.  J.  Num. 
Meth.  Engrg.,  27,  pp.  361-382,  (1989). 

[21]  Tan,  S.  C.,  “A  Progressive  Failure  Model  of  Composite  Laminates  Containing  Openings,”  J. 
Comp.  Mats.  25,  pp.  556-577,  (1991). 

[22]  Ochoa,  0.  0.  and  Engblom,  J.  J.,  “Analysis  of  Progressive  Failure  in  Composites,”  Comp. 
Sci.  Tech.,  28,  pp.  87-102,  (1987). 

[23]  Lee,  J.  D.,  “Three  Dimensional  Finite  Element  Analysis  of  Damage  Accumulation  in  Com¬ 
posite  Laminate,”  Comp.  Struct,  Vol.  15,  No.  3,  pp.  335-350,  (1982). 

[24]  Sun,  C.  T.,  “Failure  Analysis  of  Laminated  Composites  by  Using  Iterative  Three  Dimensional 
Finite  Element  Method,”  Comp.  Struct,  Vol.  33,  No.  1,  pp.  41-47,  (1989). 

[25]  Carriveau,  G.  W.,  “Benchmarking  of  the  State-of-the-Art  in  Nondestructive  Test¬ 
ing/Evaluation  for  Applicability  in  the  Composite  Armored  Vehicle  (CAV)  Advanced  Tech¬ 
nology  Demonstrator  (ATD)  program,”  TARDEC  TC  No.  13604,  Nov.  1993. 

[26]  Adams,  R.  D.  and  Cawley,  P.,  “A  Review  of  Defect  Types  and  Nondestructive  Testing  Tech¬ 
niques  for  Composites  and  Bonded  Joints,”  NDT  International,  Vol.  21,  No.  4,  p.  208-222, 
(1988). 

[27]  Rogovsky,  A.  J.,  “Engineering  Selection  of  NDT  Techniques  and  Accept/Reject  Criteria  for 
Composites,”  New  Direction  in  the  Nondestructive  Evaluation  of  Advanced  Materials.,  pp. 
45-52,  (1988). 

[28]  Murri,  W.  J.,  Sermon,  B.  W.,  Andersen,  R.  N.,  Martinez,  L.  A.  and  VanderHeiden,  E.  J., 
“Defects  in  Thick  Composites  and  some  Methods  to  Locate  Them,”  Review  of  Progress  in 
Quantitative  Nondestructive  Evaluation.,  10(B),  pp.  1583-1590,  (1991). 

[29]  Berthelot,  J.  M.  and  Rhazi,  J.,  “Acoustic  Emission  in  Carbon  Fibre  Composites,”  Comp.  Sci. 
Tech.,  37,  pp.  411-428,  (1990). 

[30]  Daniel,  I.  M.  and  Wooh,  S.  C.,  “Ultrasonic  Characterization  of  Defects  and  Damage  in  Thick 
Composites,”  Review  of  Progress  in  Quantitative  Nondestructive  Evaluation.,  9(B),  pp.  1489- 
1496, (1990). 


43 


[31]  Hamstad,  M.  A.,  “A  Review:  Acoustic  Emission,  a  Tool  for  Composite  Material  Studies,” 
Exp.  Mech.,  pp.  7-13,  (1986). 

[32]  Michaels,  T.  E.  and  Davidson,  B.  D.,  “Ultrasonic  Inspection  Detects  Hidden  Damage  in 
Composites,”  Advanced  Materials  &  Processes.,  143(3),  pp.  34-38,  (1993). 

[33]  Hsu,  D.  K.  and  Margetan,  F.  J.,  “Analysis  of  Acousto- Ultrasonic  Signal  in  Unidirectional 
Thick  Composites  Using  the  Slowness  Surfaces,”  J.  Comp.  Mats.,  26(4),  pp.  1050-1061, 
(1992). 

[34]  Schechter,  R.  S.,  Chaskelis,  H.  H.,  Mignogna,  R.  B.  and  Delsanto,  P.  P.,  “Real-Time  Parallel 
Computation  and  Visualization  of  Ultrasonic  Pulses  in  Solids,”  Science,  Vol.  265,  pp.  1188- 
1192, (1994). 

[35]  Bakuckas,  J.  G.,  Prosser,  W.  H.  and  Johnson,  W.  S.,  “Monitoring  Damage  Growth  in  Tita¬ 
nium  Matrix  Composites  Using  Acoustic  Emission,”  J.  Comp.  Mats.,  28,  pp.  305-316,  (1994). 

[36]  Hobbs,  C,  and  Temple,  A.,  “The  Inspection  of  Aerospace  Structures  Using  Transient  Ther¬ 
mography,”  British  Journal  of  NOT.,  35(4),  pp.  183-189,  (1993). 

[37]  Jones,  T.  and  Berger,  H.,  “Thermographic  Detection  of  Impact  Damage  in  Graphite-Epoxy 
Composites,”  Mater.  Eval.,  50(12),  pp.  1446-1453,  (1992). 

[38]  Henry,  D.  L.,  McCall,  C.  D.  and  Burns,  B.  P.,  “Nondestructive  Evaluation  of  Composite  Ma¬ 
terials  by  Computed  Tomography  (CT),”  Army  Research  Laboratory,  ARL-MR-32,  (1993). 

[39]  Lambrineas,  P.,  Davis,  J.  R.,  Suendermann,  B.,  Wells,  P.,  Thomson,  K.  R.,  Woodward,  R.  L., 
Egglestone,  G.  T.  and  Challis,  K.,  “X-Ray  Computed  Tomography  Examination  of  Inshore 
Mine-Hunter  HuU  Composite  Material,”  NDT  &  E  International,  pp.  207-213,  (1991). 

[40]  Doering,  E.  R.,  Basart,  J.  P.  and  Bray,  J.  N.,  “Three-Dimensional  Flaw  Reconstruction  and 

Dimensional  Analysis  Using  a  Real-Time  X-Ray  Imaging  System,”  NDT  &  E  International, 

Vol.  26,  No.  1,  pp.  7-17,  (1993). 

✓ 

[41]  Frock,  B.  G.  and  Martin,  R.  W.,  “Applications  of  Digital  Image  Enhancement  Techniques 
for  Improved  Ultrasonic  Imaging  of  Defects  in  Composite  Materials,”  Review  of  Progress  in 
Quantitative  Nondestructive  Evaluation.,  6(A),  pp.  781-789,  (1986). 

[42]  Potet,  P.,  Jeannin,  P.  and  Bathias,  C.,  “The  Use  of  Digital  Image  Processing  in  Vibrother- 
mographic  Detection  of  Impact  Damage  in  Composite  Materials,”  Mater.  Eval.,  45(4),  pp. 
466-470,  (1987). 

[43]  Steiner,  K.  V.,  “Image  Enhancement  Techniques  for  Ultrasonic  NDE  Applications,”  Compos¬ 
ite  Materials:  Testing  and  Design  (10th  Volume).,  ASTM  STP  1120,  pp.  330-343,  (1992). 

[44]  Webster,  R.  T.,  Das,  P.  and  Werner,  R.,  “Ultrasonic  Imaging  for  Nondestructive  Evaluation  of 
Composite  Materials  with  Digital  Image  Enhancement,”  1980  Ultransonic  Symposium  Proc., 
VoL  2,  pp.  873-876,  (1980). 

[45]  Weight,  K.,  “An  Overview  on  NDE  Methods  for  Thick  Composites  and  a  Proposal  for  Analysis 
of  Computed  Tomography  Data,”  U.S.  Army  Research  Laboratory,  ARL  TR  516,  (1994). 

[46]  Chou,  S.  C.  and  DeLuca,  E.,  ”Dynamic  Response  of  S-2  Glass  Reinforced  Plastic  Structural 
Armor,”  U.S.  Army  Research  Laboratory,  ARL-SR-5,  (1993). 


44 


[47]  Hughes,  J.  D.  H.,  “The  Carbon  Fibre/Epoxy  Interface  -  a  Review,”  Comp.  Sci.  Tech.,  Vol. 
41,  pp.  13-45,  (1991). 

[48]  TsangaraJcis,  N.  and  Taleghani,  B.  K.,  “Biaxial  Flexural  Testing  of  a  Carbon  Fiber  Reinforced 
Epoxy  Composite,”  U.S.  Army  Research  Laboratory,  ARL-TR-221,  (1993). 

[49]  Dreumel,  W.  H.  M.  and  Kamp,  J.  L.  M.,  “Non  Hookean  Behaviour  in  the  Fibre  Direction  of 
Carbon  Fibre  Composites  and  the  Influence  of  Fibre  Waviness  on  the  Tensile  Properties,”  J. 
Comp.  Mats.,  11,  pp.  461-469,  (1977). 

[50]  Hahn,  H.  T.  and  Tsai,  S.  W.,  “Nonlinear  Elastic  Behavior  of  Unidirectional  Composite  Lam¬ 
inate,”  J.  Comp.  Mats.,  7,  pp.  102-110,  (1973). 

[51]  Amijima,  S.  and  Adachi,  T.,  “Nonlinear  Stress-Strain  Response  of  Laminated  Composites,” 
J.  Comp.  Mats.,  13,  pp.  206-218,  (1979). 

[52]  Sims,  D.  F.,  “In-Plane  Shear  Strain  Response  of  Unidirectional  Composite  Materials,”  J. 
Comp.  Mats.,  7,  pp.  124-130,  (1973). 

[53]  Hashin,  Z.,  “Failure  Criteria  for  Unidirectional  Fiber  Composites,”  J.  Appl.  Mech.,  47,  pp. 
329-334,  (1980). 

[54]  Collins,  J.  A.,  Failure  of  Materials  in  Mechanical  Design,  John  Wiley  &  Sons,  (1981). 

[55]  Hill,  R.,  ‘A  Theory  of  Yielding  and  Plastic  Flow  of  Anisotropic  Metals,”  Proceedings  of  the 
Royal  Society,  Series  A,  193,  p.  281,  (1948). 

[56]  Tsai,  S.  W.,  “Strength  Characterization  of  Composite  Materials,”  NASA  CR-224,  (1965). 

[57]  Hoffman,  0.,  “The  Brittle  Strength  of  Orthotropic  Materials,”  J.  Comp.  Mats.,  1,  p.  200, 
(1967). 

[58]  Tsai,  S.  W.  and  Wu,  E.  M.,  “A  General  Theory  of  Strength  for  Anisotropic  Materials,”  J. 
Comp.  Mats.,  5,  p.  58,  (1971). 

[59]  Sandhu,  R.  S.,  “A  Survey  of  Failure  Theories  of  Isotropic  and  Anisotropic  Materials,”  Tech¬ 
nical  Report,  AFFDL-TR- 72-71. 

[60]  Tsai,  S.  W.,  “A  Survey  of  Macroscopic  Failure  Criteria  for  Composite  Materials,”  Technical 
Report,  AFWAL-TR-84-4025. 

[61]  Reddy,  N.  J.  and  Pandey,  A.  K.,  “A  First-Ply  Failure  Analysis  of  Composite  Laminates,” 
Comp.  Struct.,  Vol.  25,  No.  3,  pp.  371-393,  (1987). 

[62]  Narayanaswami,  R.,  “Evaluation  of  the  Tensor  Plynomial  and  Hoffman  Strength  Theories  for 
Composite  Materials,”,  J.  Comp.  Mats.,  11,  pp.  366-377,  (1977). 

[63]  Wu,  R.  Y.  and  Stachursky,  Z.,  “Evaluation  of  the  Normal  Stress  Interaction  Parameter  in 
the  Tensor  Polynomial  Strength  Theory  for  Anisotropic  Materials,”  J.  Comp.  Mats.,  18,  pp. 
456-463,  (1984). 

[64]  Wu,  R.  Y,  and  Stachursky,  Z.,  “Evaluation  of  the  Normal  Stress  Interaction  Parameter  in 
the  Tensor  Polynomial  Strength  Theory  for  Anisotropic  Materials,”  J.  Comp.  Mats.,  18,  pp. 
456-463, (1984). 


45 


[65]  Chang,’  F.  K.  and  Lessard,  L.  B.,  “Damage  Tolerance  of  Laminated  Composites  Containing 
an  Open  Hole  and  Subjected  to  Compressive  Loading:  Part  I  -  Analysis,”  J.  Comp.  Mats. ,2^, 
pp.  2-43,  (1991). 

[66]  Choi,  H.  Y.,  Wu,  H.  T.  and  Chang,  F.,  “A  New  Approach  Toward  Understanding  Damage 
Mechanisms  and  Mechanics  of  Laminated  Composites  Due  to  Low  Velocity  impact:  Part  II 
-  Analysis,”  J.  Comp.  Mats.,  Vol.  25,  pp.  1012-1038,  (1991). 

[67]  Lagoudas,  D.  C.,  and  Saleh,  A.  M.,  “Compressive  Failure  Due  to  Kinking  of  Fibrous  Com¬ 
posites”,  J.  Comp.  Mats.,  27.  pp.  83-106,  (1993). 

[68]  Budiansky,  B.,  “Micromechanics”,  Comp.  Struct.,  Vol.  16,  No.  1-4,  pp.  3-12,  (1983). 

[69]  Gao,  Z.  and  Reifsnider,  K.  L.,  “Tensile  Failure  of  Composites:  Influence  of  Interface  and 
Matrix  Yielding,”  J.  Comp.  Tech.  Res.,  Vol.  14,  No.  4,  pp.  201-210,  (1992). 

[70]  Chang,  F.  K.  and  Chen,  M.,  “The  In  Situ  Ply  Shear  Strength  Distribution  in  Graphite/Epoxy 
Laminated  Composites,”  J.  Comp.  Mats.,  21,  pp.  708-732,  (1987). 

[71]  Yamada,  S.  E.  and  Sun,  C.  T.,  “Analysis  of  Laminate  Strength  and  its  Distribution,”  J. 
Comp.  Mats.,  12,  pp.  275-284,  (1978). 

[72]  Chang,  F.  K.,  Scott,  R.  A.  and  Springer,  G.  S.,  “The  Effect  of  Laminate  Cbnflguration  on 
Characteristic  Lengths  and  Rail  Shear  Strength,”  J.  Comp.  Mats.,  18,  pp.  290-296,  (1984). 

[73]  Flaggs,  D.  L.  and  Kural,  M.  H.,  “Experimental  Determination  of  In  Situ  Transverse  Lamina 
Strength  in  Graphite/Epoxy  Laminate,”  J.  Comp.  Mats.,  16,  pp.  103-116,  (1982). 

[74]  Peter,  P.  W.  M.,  “The  Strength  Distribution  of  90°  Plies  in  0/90/0  Graphite  Epoxy  Lami¬ 
nates,”  J.  Comp.  Mats.,  18,  pp.  545-556,  (1984). 

[75]  Garrett,  K.  W.  and  Bailey,  J.  E.,  “Multiple  Transverse  Fracture  in  90°  Cross-Ply  Laminates 
of  a  Glass  Fibre-Reinforced  Polyester,”  J.  Mater.  Sci.,  pp.  157-168,  (1977). 

[76]  Talreja,  R.,  “Transverse  Cracking  and  Stiffness  Reduction  in  Composite  Laminates,”  J.  Comp. 
Mats.,  19,  pp.  355-375,  1985. 

[77]  Highsmith,  A.  L.  and  Reifsnider,  K.  L.,  “Stiffness- Reduction  Mechanisms  in  Composite  Lam¬ 
inates,”  in  Damage  in  Composite  Materials,  ASTM  STP  775,  American  Society  for  Testing 
and  Materials,  Philadelphia,  PA,  pp.  103-117,  (1982). 

[78]  Kistner,  M.  D.,  Whitney,  J.  M.  and  Browning,  C.  E.,  “First  Ply  Failure  of  Graphite/Epoxy 
Laminates,”  Poc.  U.S.-Japan  Conference  on  Composite  Materials,  NASA-Langley,  Hampton, 
Virginia,  (June  1983). 

[79]  Doucet,  A.  B.  and  Qin,  J.,  “Microstructural  Analysis  of  Impact  Damage  in  Two  Glass/Epoxy 
Composite  Laminates,”  Comp.  Engrg.,  3,  pp.  55-68,  (1993). 

[80]  Rosen,  B.  W.,  “Mechanics  of  Composite  Strengthening,”  Fiber  Composite  Materials,  Chapter 
3,  American  Society  for  Metals,  (1965). 

[81]  Tsai,  S.  W.  and  Hahn,  H.  T.,  Introduction  to  Composite  Materials,  Lancaster  PA:  Technomic 
Publishing  Company,  (1980). 


46 


[82]  Rosen,  B.  W.,  “Tensile  Failure  of  Fibrous  Composites,”  AIAA,  11,  pp.  1985-1991,  (1964). 

[83]  Nuismer,  R.  J.  and  Tan,  S.  C.,  “Constitutive  Relations  of  a  Cracked  Composite  Lamina,”  J. 
Comp.  Mats.,  22,  pp.  306-321,  (1988). 

[84]  Pyrz,  R.,  “A  Micromechanically-Based  Model  for  Composite  Materials  with  Matrix  Cracks,” 
Comp.  Engrg.,  Vol.  2,  No.  8,  pp.  619-629,  (1992). 

[85]  Jones,  R.  M.,  Mechanics  of  Composite  Materials,  McGRAW-HILL,  (1975). 

J  [86]  O’Brien,  T.  K.,  “Characterization  of  Delamination  Onset  and  Growth  in  a  Composite  Lami¬ 

nate,”  Damage  in  Composite  Materials,  ASTM  STP  775,  American  Society  for  Testing  and 
Materials,  pp.  140-167,  (1982). 

[87]  Kim,  R.  Y.  and  Soni,  S.  R.,  “Experimental  and  Analytical  Studies  on  the  Onset  of  Delami¬ 
nation  in  Laminated  Composites,”  J.  Comp.  Mats.,  18,  pp.  70-80,  (1984). 

[88]  Kim,  R.  Y.  and  Soni,  S.  R.,  “Failure  of  Composite  Laminates  due  to  Combined  Interlaminar 
Shear,”  Composites  ‘86:  Recent  Advances  in  Japan  and  the  United  States,  pp.  341-350,  (1986). 

[89]  Brewer,  J.  C.  and  Lagace,  P.  A.,  “Quadratic  Stress  Criterion  for  Initiation  of  Delamination,” 
J.  Comp.  Mats.,  22,  pp.  1141-1155,  (1988). 

[90]  Lui,  D.,  “Impact-Induced  Delaminations  -  A  View  of  Bending  Stiffness  Mismatching,”  J. 
Comp.  Mats.,  24,  pp.  674-692,  (1988). 

[91]  Nairn,  J.  A.,  Lui,  S.  and  Chen,  H.,  “Longitudinal  Splitting  in  Epoxy  and  K-Polymer  Com¬ 
posites:  Shear  Lag  Analysis  Including  the  Effect  of  Fiber  Bridging,”  J.  Comp.  Mats.,  25,  pp. 
1086-1107,  (1991). 

[92]  Dym,  C.  L.  and  Shames,  I.  H.,  Solid  Mechanics-A  Variational  Approach,  McGraw  Hill,  New 
York,  (1973). 

[93]  Rybicki,  E.  F,.  Schmueser,  D.  W.  and  Fox,  J.,  “An  Energy  Release  Rate  Approzich  for  Stable 
Crack  Growth  in  the  Free-Edge  Delamination  Problem,”  J.  Comp.  Mats.,  11,  pp.  470-487, 
(1977). 

[94]  Wang,  S.  S.,  “Fracture  Mechanics  for  Delamination  Problems  in  Composite  Materials,”  J. 
Comp.  Mats.,  17,  pp.  210-223,  (1983). 

>  [95]  Rybicki,  E.  F.  and  Kanninen,  M.  F.,  “A  Finite  Element  Calculation  of  Stress  Intensity  Factors 

by  a  Modified  Crack  Closure  Integral,”  Engrg.  Fract.  Mech.,  Vol.  9,  No.  4,  pp.  931-938,  (1977). 

,,  [96]  Irwin,  G.  R.,  Handbuch  der  Physik,  6,  p.  551,  (1958). 

[97]  0.  C.  Zienkiewicz  and  R.  L.  Taylor,  The  Finite  Element  Method,  Vol.  I,  Basic  Introduction 
and  Linear  Problem,  4th  edn.  McGraw-Hill  (1989). 

[98]  Mindlin,  R.  D.,  “Influence  of  Rotary  Inertia  and  Shear  on  Flexural  Motions  of  Isotropic, 
Elastic  Plates,”  J.  Appl.  Mech.  ASME,  18,  pp.  31-38,  (1951). 

[99]  Lo,  K.  H.,  Christiansen,  R.  M.  and  Wu,  E.  M.,  “A  Higher-Order  theory  of  Plate  Deformation, 
Part  1:  Homogeneous  Plates  and  Part  2:  Laminated  Plates,”  J.  Appl.  Mech.  ASME,  44,  pp. 
663-676, (1977). 


47 


[100]  Reissner,  E.,  “Reflections  on  the  Theory  of  Elastic  Plates,”  Appl.  Mech.  Rev.,  38,  pp.  1453- 
1464, (1985). 

[101]  Reddy,  J.  N.,  “On  Refined  Computational  Models  of  Composite  Laminates,”  Int.  J.  Numer. 

Meth.  Engr.,  27,  pp.  361-382,  (1989). 

[102]  Noor,  A.  K.  and  Burton,  W.  S.,  “Assessment  of  Shear  Deformable  Theories  for  Multilayered 
Composite  Plates,”  Appl.  Mech.  Rev.,  42,  pp.  1-12,  (1989). 

[103]  Tessler,  A.  and  Saether,  E.,  “A  Computationally  Viable  Higher-Order  Theory  for  Laminated 

Composite  Plates,”  Int.  J.  Num.  Meth.  Engr.,  31,  pp.  1069-1086,  (1991).  ‘ 

[104]  Sciuva,  M.  D.,  “Bending,  Vibration  and  Buckling  of  Simply  Supported  Thick  Multilayered 

Orthotropic  Plates:  An  Evaluation  of  a  New  Displacement  Model,”  J.  Sound  and  Vibration,  , 

105,  pp.  425-442,  (1986). 

[105]  Robbins,  D.  H.  and  Reddy,  N.  J.,  “Modelling  of  Thick  Composites  Using  a  Layerwise  Lami¬ 
nate  Theory,”  Int.  J.  Numer.  Meth.  Engr.,  36,  pp.  655-677,  (1993). 

[106]  T.  H.  H.  Pian  and  P.  Tong,  “Basis  of  Finite  Element  Methods  for  Solid  Continua,”  Int.  J. 

Num.  Meth.  Engrg.  1,  pp.  3-28,  (1969). 

[107]  T.  H.  H.  Pian  and  D.  P.  Chen,  “Alternative  Ways  for  Formulation  of  Hybrid  Stress  Elements,” 

Int.  J.  Num.  Meth.  Engrg.  18,  pp.  1679-1684,  (1982). 

[108]  T.  H.  H.  Pian  and  D.  P.  Chen,  “On  the  Suppression  of  Zero  Energy  Deformation  Modes,” 

Int.  J.  Num.  Meth.  Engrg.  19,  pp.  1741-1752  (1983). 

[109]  T.  H.  H.  Pian  and  K.  Sumihara,  “Rational  Approach  for  Assumed  Stress  Finite  Element,” 

Int.  J.  Num.  Meth.  Engrg.  20,  pp.  1685-1695,  (1984). 

[110]  T.  H.  H.  Pian  and  P.  Tong,  “Relations  Between  Incompatible  Displacement  Model  and  Hybrid 
Stress  Model,”  Int.  J.  Num.  Meth.  Engrg.  22,  pp.  173-181,  (1986). 

[111]  K.  Y.  Sze  and  C.  L.  Chow,  “Efficient  Hybrid/Mixed  Elements  Using  Admissible  Matrix 
Formulation,”  Comp.  Meth.  Applied  Mech.  and  Engrg.,  99,  pp.  1-26  (1992). 

[112]  K.  Y.  Sze,  “Efficient  Formulation  of  Robust  Hybrid  Elements  Using  Orthogonal  Stress/Strain 
Interpolants  and  Admissible  Matrix  Formulation,”  Int.  J.  Num.  Meth.  Engrg.  35,  pp.  1-20, 

(1992).  ^ 

[113]  Bathe,  K.  J.,  Finite  Element  Procedures  in  Engineering  Analysis,  Prentice-Hall,  (1982). 

[114]  C.  Lanczos,  Applied  Analysis,  Dover  Publications,  (1988).  -f 

[115]  Cook,  R.  D,  Concepts  and  Applications  of  Finite  Element  Analysis,  2'^'^ Edition,  John  Wiley 
&  Sons,  (1981). 

[116]  Press,  W.  H.,  Teukolsky,  S.  A.,  Vetterlin,  W.  T.  and  Flannery,  B.  P.,  Numerical  Recipes, 

2^^^ Edition,  Cambridge  University  Press,  (1992). 

[117]  Larsson,  P.  L.,  “On  Multiple  Delamination  Buckling  and  Growth  in  Composite  Plates,”  Int. 

J.  Solids  Structures,  Vol.  27,  No.  13,  pp.  1623-1637,  (1991). 


48 


[118]  Staxakers,  B.  and  Andersson,  B.,  “Nonlinear  Plate  Theory  Applied  to  Delamination  in  Com¬ 
posites,”  J.  Mech.  Phys.  Solids',  36,  pp.  689-718,  (1988). 

[119]  Jeusette,  J.  P.  and  Sonzogni,  V.,  “A  Projected  Conjugate  Gradient  Method  for  Structural 
Stability  Analysis  with  Linear  Constraints,”  Comp.  Struct,,  Vol.  33,  No.  1,  pp.  31-39,  (1989). 

[120]  Whitcomb,  J.  D.,  “Analysis  of  a  Laminate  with  a  Postbuckled  Embedded  Delamination, 
Including  Contact  Effects,”  J.  Comp.  Mats.,  26,  pp.  1523-1535,  (1992). 


7 


49 


No.  of 
Copies 


DISTRIBUTION  LIST 


To 


1  Office  of  the  Under  Secretary  of  Defense  for  Research  and  Engineering,  The  Pentagon, 
Washington,  DC  20301 

Director,  U.S.  Army  Research  Laboratory,  2800  Powder  Mill  Road,  Adelphi,  MD  20783-1197 
1  ATTN:  AMSRL-OP-SD-TP,  Technical  Publishing  Branch 
1  AMSRL-OP-SD-TA,  Records  Management 

1  AMSRL-OP-SD-TL,  Technical  Library 

Commander,  Defense  Technical  Information  Center,  Cameron  Station,  Building  5, 

5010  Duke  Street,  Alexandria,  VA  23304-6145 

2  ATTN:  DTIC-FDAC 

1  MIA/CINDAS,  Purdue  University,  2595  Yeager  Road,  West  Lafayette,  IN  47905 

Commander,  Army  Research  Office,  P.O.  Box  12211,  Research  Triangle  Park, 

NC  27709-2211 

1  ATTN:  Information  Processing  Office 

Commander,  U.S.  Army  Materiel  Command,  5001  Eisenhower  Avenue,  Alexandria,  VA  22333 
1  ATTN:  AMCSCI 


Commander,  U.S.  Army  Materiei  Systems  Analysis  Activity,  Aberdeen  Proving  Ground, 
MD  21005 

1  ATTN:  AMXSY-MP,  H.  Cohen 

Commander,  U.S.  Army  Missile  Command,  Redstone  Arsenal,  AL  35809 
1  ATTN:  AMSMI-RD-CS-R/Doc 

Commander,  U.S.  Army  Armament,  Munitions  and  Chemical  Command,  Dover,  NJ  07801 
1  ATTN:  Technical  Library 

Commander,  U.S.  Army  Natick  Research,  Development  and  Engineering  Center 
Natick,  MA  01760-5010 
1  ATTN:  SATNC-MI,  Technical  Library 

Commander,  U.S.  Army  Satellite  Communications  Agency,  Fort  Monmouth,  NJ  07703 
1  ATTN:  Technical  Document  Center 

Commander,  U.S.  Army  Tank-Automotive  Command,  Warren,  Ml  48397-5000 
1  ATTN:  AMSTA-ZSK 
1  AMSTA-TSL,  Technical  Library 

President,  Airborne,  Electronics  and  Special  Warfare  Board,  Fort  Bragg,  NC  28307 
1  ATTN:  Library 

Director,  U.S.  Army  Research  Laboratory,  Weapons  Technology,  Aberdeen  Proving  Ground 
MD  21005-5066 
1  ATTN:  AMSRL-WT 


50 


No.  of 
Copies 


To 


Commander,  Dugway  Proving  Ground,  UT  84022 
1  ATTN:  Technical  Library,  Technical  Information  Division 

Commander,  U.S.  Army  Research  Laboratory,  2800  Powder  Mill  Road,  Adelphi,  MD  20783 
1  ATTN:  AMSRL-SS 

Director,  Benet  Weapons  Laboratory,  LCWSL,  USA  AMCCOM,  Watervliet,  NY  12189 
1  ATTN:  AMSMC-LCB-TL 

1  AMSMC-LCB-R 

1  AMSMC-LCB-RM 

1  AMSMC-LCB-RP 

Commander,  U.S.  Army  Foreign  Science  and  Technology  Center,  220  7th  Street,  N.E., 
Charlottesville,  VA  22901-5396 

3  ATTN:  AIFRTC,  Applied  Technologies  Branch,  Gerald  Schlesinger 

Commander,  U.S.  Army  Aeromedical  Research  Unit,  P.O.  Box  577,  Fort  Rucker,  AL  36360 
1  ATTN:  Technical  Library 

U.S.  Army  Aviation  Training  Library,  Fort  Rucker,  AL  36360 
1  ATTN:  Building  5906-5907 

Commander,  U.S.  Army  Agency  for  Aviation  Safety,  Fort  Rucker,  AL  3636 
1  ATTN:  Technical  Library 

Commander,  Clarke  Engineer  School  Library,  3202  Nebraska  Ave.,  N.,  Fort  Leonard  Wood, 
MO  65473-5000 
1  ATTN:  Library 

Commander,  U.S.  Army  Engineer  Waterways  Experiment  Station,  P.O.  Box  631,  Vicksburg, 
MS  39180 

1  ATTN:  Research  Center  Library 

Commandant,  U.S.  Army  Quartermaster  School,  Fort  Lee,  VA  23801 
1  ATTN:  Quartermaster  School  Library 

Naval  Research  Laboratory,  Washington,  DC  20375 
1  ATTN:  Code  6384 

Chief  of  Naval  Research,  Arlington,  VA  22217 
1  ATTN:  Code  471 

Commander,  U.S.  Air  Force  Wright  Research  and  Development  Center,  Wright-Patterson 
Air  Force  Base,  OH  45433-6523 
1  ATTN:  WRDC/MLLP,  M.  Forney,  Jr. 

1  WRDC/MLBC,  Mr.  Stanley  Schulman 

U.S.  Department  of  Commerce,  National  Institute  of  Standards  and  Technology,  Gaithersbura 
MD  20899 

1  ATTN:  Stephen  M  Hsu,  Chief,  Ceramics  Division,  Institute  for  Materials  Science 
and  Engineering 


51 


Committee  on  Marine  Structures,  Marine  Board,  National  Research  Council,  2101  Constitution 
Avenue,  N.W.,  Washington,  DC  2041  8 

Materials  Sciences  Corporation,  Suite  250,  500  Office  Center  Drive,  Fort  Washington, 

PA  19034 

Charles  Stark  Draper  Laboratory,  555  Technology  Square,  Cambridge,  MA  02139 

General  Dynamics,  Convair  Aerospace  Division,  P.O.  Box  748,  Fort  Worth,  TX  76101 
ATTN:  Mfg.  Engineering  Technical  Library 

Plastics  Technical  Evaluation  Center,  PLASTEC,  ARDEC,  Bldg.  355N,  Picatinny  Arsenal, 

NJ  07806-5000 
ATTN:  Harry  Pebly 

Department  of  the  Army,  Aerostructures  Directorate,  MS-266,  U.S.  Army  Aviation  R&T 
Activity  -  AVSCOM,  Langley  Research  Center,  Hampton,  VA  23665-5225 

NASA  -  Langley  Research  Center,  Hampton,  VA  23665-5255 

U.S.  Army  Vehicle  Propulsion  Directorate,  NASA  Lewis  Research  Center, 

2100  Brookpark  Road,  Cleveland,  OH  44135-3191 
ATTN:  AMSRL-VP 

Director,  Defense  Intelligence  Agency,  Washington,  DC  20340-6053 
ATTN:  0DT-5A,  Mr.  Frank  Jaeger 

U.S.  Army  Communications  and  Electronics  Command,  Fort  Monmouth,  NJ  07703 
ATTN:  Technical  Library 

U.S.  Army  Research  Laboratory,  Electronic  Power  Sources  Directorate, 

Fort  Monmouth,  NJ  07703 
ATTN:  Technical  Library 

Director,  U.S.  Army  Research  Laboratory,  Watertown,  MA  02172-0001 
ATTN:  AMSRL-OP-WT-IS,  Technical  Library 
Author 


